# Heart rate variability helps classify phenotype in systemic sclerosis

Stéphane Delliaux, Abdou Khadir Sow, Anass Echcherki, Audrey Benyamine, Quentin Gomes de Pinho, Fabienne Brégeon, Brigitte Granel

► To cite this version:

Stéphane Delliaux, Abdou Khadir Sow, Anass Echcherki, Audrey Benyamine, Quentin Gomes de Pinho, et al.. Heart rate variability helps classify phenotype in systemic sclerosis. Scientific Reports, 2024, 14 (1), pp.11151. 10.1038/s41598-024-60553-1 . hal-04927426

HAL Id: hal-04927426

https://hal.science/hal-04927426v1

Submitted on 17 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for updates

OPEN

# Heart rate variability helps classify phenotype in systemic sclerosis

Stéphane Delliaux<fnref n="1" />,²,³,⁴✉, Abdou Khadir Sow²,⁵, Anass Echcherki⁴, Audrey Benyamine<fnref n="1" />,⁶, Quentin Gomes de Pinho<fnref n="1" />,⁶, Fabienne Brégeon²,⁷ & Brigitte Granel<fnref n="1" />,⁶

We aimed to develop a systemic sclerosis (SSc) subtypes classifier tool to be used at the patient's bedside. We compared the heart rate variability (HRV) at rest (5-min) and in response to orthostatism (5-min) of patients (n = 58) having diffuse (n = 16, dcSSc) and limited (n = 38, lcSSc) cutaneous forms. The HRV was evaluated from the beat-to-beat RR intervals in time-, frequency-, and nonlinear-domains. The dcSSc group differed from the lcSSc group mainly by a higher heart rate (HR) and a lower HRV, in decubitus and orthostatism conditions. Stand-up maneuver lowered HR standard deviation (sd_HR), the major axis length of the fitted ellipse of Poincaré plot of RR intervals (SD2), and the correlation dimension (CorDim) in the dcSSc group while increased these HRV indexes in the lcSSc group (p = 0.004, p = 0.002, and p = 0.004, respectively). We identified the 5 most informative and discriminant HRV variables. We then compared 341 classifying models (1 to 5 variables combinations × 11 classifier algorithms) according to mean squared error, logloss, sensitivity, specificity, precision, accuracy, area under curve of the ROC-curves and F1-score. F1-score ranged from 0.823 for the best 1-variable model to a maximum of 0.947 for the 4-variables best model. Most specific and precise models included sd_HR, SD2, and CorDim. In conclusion, we provided high performance classifying models able to distinguish diffuse from limited cutaneous SSc subtypes easy to perform at the bedside from ECG recording. Models were based on 1 to 5 HRV indexes used as nonlinear markers of autonomic integrated influences on cardiac activity.

**Keywords** Systemic sclerosis, Heart rate variability, Nonlinear dynamics, Cardiovascular function, Machine learning

Systemic sclerosis (SSc) is a rare autoimmune disease combining vascular abnormalities with skin and deep organs fibrotic damaging that leads to the highest mortality among rheumatic diseases.

SSc is a heterogeneous disease and many efforts have been made for helping physicians to determinate patients' prognosis. Since 1988, a sub-classification based on sclerosis skin extension assessed by palpation of 17 areas allows the diagnosis of the diffuse cutaneous SSc (dcSSc) and limited cutaneous SSc (lcSSc). Natural history, rate of progression, severity of organ damage and dysfunction, and mortality rate of disease depend in part on these 2 SSc subtypes. The dcSSc having the worst prognosis, identifying as soon as possible and with a robust manner SSc subtype remains critical for risk stratification and SSc patients care strategy<fnref n="1" />.

Clinical and biological evaluation has largely driven SSc risk stratification over the past decades, however, integration of new markers might improve disease phenotyping. As asymptomatic but also severe cardiac abnormalities, including myocarditis, congestive heart failure, arrhythmia, focal fibrosis, and impaired ventricular relaxation, are associated with systemic sclerosis and play a significant role in a considerable proportion of related fatalities<fnref n="1" />, markers of cardiac activity have been studied. Heart rate variability (HRV), that is often considered to reflect the autonomic nervous system (ANS) modulation of the sinus node rhythmic activity, has been proposed to assess cardiovascular function impairment in SSc since 1994². Several studies have reported a lower spontaneous HRV in SSc patients compared to healthy controls in supine resting conditions³,⁴. The HRV index that is low-frequency spectral power has been associated with an enhanced resting sympathetic hyperactivity⁵. Additionally, SSc patients experienced first, a reduced heart rate (HR) response to tilt test highlighting impaired baroreflex modulation of the cardiovascular autonomic control and second a correlated high-frequency spectral power—lung function suggesting a respiratory sinus arrhythmia-supported HRV⁵,⁶. But HRV is also and more

<footnote n="1">¹INSERM, INRAE, C2VN, Aix Marseille Univ, Marseille, France. ²Explorations Fonctionnelles Respiratoires, AP-HM, Hôpital Nord, Marseille, France. ³CNRS, CPT, Aix Marseille Univ, Marseille, France. ⁴Laënnec Institute – Digital Sciences for Health, Aix Marseille Univ, Marseille, France. ⁵Laboratoire de Physiologie, Cheikh Anta Diop University, Dakar, Senegal. ⁶Service de Médecine Interne, AP-HM, Hôpital Nord, Marseille, France. ⁷AP-HM, Microbes Evolution Phylogeny and Infections (MEPHI), IHU-Méditerranée Infection, Aix Marseille Univ, Marseille, France. ✉email: stephane.delliaux@univ-amu.fr</footnote>
newly considered as the result of converging common and interacting sinus activity drivers that are largely interactive⁷. These drivers include, but non-exclusively, sympathetic and parasympathetic spontaneous and reflexively-driven tones, respiration pattern through frequency and volume activating pulmonary stretch reflex, non-vagal intrathoracic mechanical effects, hypercapnia, vagal-mediated arterial baroreflex effects, and local cardiac nodal (peptides) influences. Besides, SSc is known to potentially accumulate disturbances of several of these sinus node activity determinants such as autonomic nervous system, baroreflex function, lung function, and thoracic mechanics impairments. It seems therefore highly challenging to address the complexity of the issue in SSc patients through the only time- and frequency- domains.

Indeed, HRV has also been used to characterize cardiovascular function in children with systemic and localized scleroderma and no difference has been shown⁸. In adults of different SSc subsets, a recent study failed to prove the interest of some HRV indexes, particularly from spectral analysis, to distinguish dcSSc and lcSSc subtypes⁹. Nevertheless, HRV metrics and more specifically nonlinear ones as symbolic dynamics and conditional entropy confirmed their ability to detect a predominant sympathetic modulation of sinus node rhythmic activity at rest and a blunted response to orthostatism stress especially at the SSc early stage⁹.

All these exploratory works highlight the interest of HRV analysis, particularly through nonlinear metrics, to study SSc impact on cardiovascular function. To date, none provide discriminant HRV markers to distinguish the different SSc subtypes. A bioinformatic marker for SSc subtypes classification is still lacking. The aim of this research was to propose a classifying tool including HRV indexes only, easy to compute, and to be used at the patient's bedside. To do so, we tested the hypothesis that a combination of few variables describing sinus node rhythmic activity and particularly describing the nonlinear properties of its dynamics as correlation dimension (CorDim), could allow to distinguish dcSSc from lcSSc. In accordance with our knowledge and experience previously acquired¹⁰ and because CorDim is typically used as an index of the overall complexity of the system dynamics estimated from a time series, CorDim should be specifically efficient in multiple physiological functions interactions situations such as that observed in SSc patients performing a stand-up test.

To meet this challenge, we compared in dcSSc and lcSSc patients HRV at rest and in response to orthostatism. We computed several linear and nonlinear HRV metrics from time-, frequency-, and nonlinear- domains as a step of a data engineering process. Several classifying models (metrics combinations) and classifying methods (classifier algorithms) have been compared to help the practitioner to choose the best classifying model according to the clinical needs.

# Methods

## Study design and participants

We did a retrospective observational study of 58 consecutive patients aged 18–74 years and suffering from SSc. Patients were identified from the consulting population of the internal medicine department of north hospital, Marseille, France, from 2018 to 2021. Power analysis was performed to determine sample size. To achieve a power of 80% (type II error) and a level of significance of 5% (two-sided type I error) for detecting a 50% group difference (i.e. > 2 vs. < 1 degrees of freedom in the RRI generating dynamical system) in expected mean value (2.2 ± 0.8, based on previous study¹⁰ and preliminary data) of the correlation dimension (CorDim), 16 patients (8 per group) were necessary. To ensure high models performances we included more than necessary subjects. Moreover, a posteriori calculation of Cohen's d metrics (equal to 1.16761) as well as Hedges' g metrics (equal to 1.455829) that is usually used for unbalanced groups showed that size effect was very strong in our study and that sample size = 76 (16/60) was sufficient.

Inclusion criteria were 18–80 years of age, diagnosis of SSc confirmed by an internist medical doctor according to the 2013 classification criteria for SSc¹¹. No specific criteria on the presence or absence of dysautonomia was checked at this step. Exclusion criteria were known heart disease, beta-blockers medications, and diabetes mellitus. Neither heart disease nor diabetes has been detected. 3 patients were identified as being treated by beta-blockers for systemic hypertension and were therefore not included. We included 60 patients and 2 were not analyzed because of incomplete data collection.

Informed consent was obtained from all subjects by the institution concerning care data exploitation for research purposes. The protocol of data exploitation was developed in accordance with the declaration of Helsinki and health data analysis approval has been provided on January 2020 by the Assistance Publique-Hôpitaux de Marseille Data Protection Officer and local ethics committee under the references PADS19-352 and WZTVF7.

## Procedures: clinical and paraclinical features

New or already known patients referred for SSc to the internal medicine department from 2018 to 2021 were retrospectively screened. All patients included in our study fulfilled the 2013 American college of Rheumatology/ European league against rheumatism classification criteria for SSc¹¹. According to the internist medical doctor specifications, patients had been explored through paraclinical biological, imaging, and functional tests including when necessary: autoantibodies tests, chest CT scan, echocardiography, right heart catheterism, lung functional tests, esophageal manometry, and cardiovascular dysautonomia checking through spectral heart rate variability assessment. The heart rate variability recordings were made prospectively on medical justifications for routine care while data were retrospectively re-analyzed with complementary technics for our research purposes. SSc patients were classified into dcSSc or lcSSc subtypes by the internist medical doctor expertise according to the 1988 LeRoy et al. classification¹². Accordingly, to the extent of skin sclerosis, two groups were considered, dcSSc and lcSSc groups.
## Heart rate variability study

The subjects were instructed to refrain from exercising, drinking alcohol, coffee or tea, and sleep deprivation for at least 24 h before their participation. The day of the HRV study, room experiment was quiet and its temperature kept constant. Subjects were asked to lie down in the examination bed and to relax with free breathing. Three lead electrocardiogram (ECG) was acquired (ECG100C, UIM100C, and MP150A-CE Biopac equipment, BIOPAC Systems Inc., Santa Barbara, CA, USA) during 15 min in supine resting condition and then during 10 min in orthostatic resting condition after stand-up. The data were digitized and stored on a personal computer for further analysis. Two periods of 5 min were extracted and analyzed: (1) 5 min of supine rest (Decubitus) after at least 5 min of stable ECG signal, and (2) 5 min of orthostatic rest (Orthostatism) starting immediately after stand-up phase (< 10 s) noise stabilization. Signals were conditioned and processed as follow and previously described¹⁰. Raw signals were sampled at 1000 Hz and digitized with a 24-bit analog-to-digital converter. Non-evident non-stationarity, such as very slow drifting of the mean or sudden changes of the variance, was observed after 3-order polynomial detrending. No additional filtering technique such as integration was used. Pan and Tompkins real-time QRS detection algorithm was used to automatically detect R-waves and build the RRI series. Tachograms were then firstly visually checked and corrected for artifacts (including ectopic beat management, i.e., detection, cancelation, and interpolation) when necessary. Then, all cardiovascular data analysis was performed from the 5 min beat-to-beat time series constituted from inter-beat intervals from ECG R waves, i.e., the RR intervals (RRI) that were extracted from raw signals. Discrete original RRI series were resampled by cubic spline interpolation with a 4 Hz sampling rate to generate equidistantly sampled time series x(i), i = 1, 2, ..., N. The HRV was assessed according to the consensual standards¹³, using Kubios software (Kubios HRV, 2.1, Biosignal Analysis and Medical Imaging Group, Kuopio, Finland), computing HRV metrics in Decubitus (_S1) as well as in Orthostatism (_S2). We also assessed stand-up induced HRV changes (_delta) calculating the difference between orthostatic and decubitus values for each variable.

## Time-domain analysis of HRV

Several classical metrics were used as indexes of the total variability that arises from both random and periodic sources¹⁴, including: the mean RRI (mean_RR) and mean HR (mean_HR), the standard deviation of the RRI (sd_RR) and of the HR (sd_HR), the square root of the mean squared differences between adjacent normal RRI (RMSSD), the count of successive normal beat lengths that differed more than 50 ms (NN50), and the percentage of successive normal beat lengths that differed more than 50 ms (pNN50).

## Frequency-domain analysis of HRV

Different metrics were used to characterize the periodic oscillations of the studied time series using the estimated power spectrum by fast Fourier transform Welch's periodogram technique. These metrics included the centroid frequency (expressed in Hz) and power (power, expressed in ms²) in 2 frequency bands of interest: high frequencies (HF, 0.15–0.4 Hz) and low frequencies (LF, 0.04–0.15 Hz). The total power (tot_power) and LH_power/HF_power ratio (LF_HF_power) were also computed. The HF_power and LF_power were also expressed as the percentage of tot_power (HF_power_perc and LF_power_perc, respectively) and normalized units (HF_power_nu and LF_power_nu, respectively). These variables were said¹³⁻¹⁵ to be indexes of the total variability of the heart rate (tot_power), the modulation of the sinus node activity by the parasympathetic component of the ANS (HF_power_nu), the modulation of the sinus node activity by both the parasympathetic and sympathetic components of the ANS (LF_power_nu), the influence of the temperature, and the sympathetic/parasympathetic balance (LF_HF_power).

## Nonlinear-domain analysis of HRV

A set of metrics are known to catch and highlight non-linear properties of the studied time series. We plotted the RRI of rank n + 1 as a function of the RRI of rank n (lag 1 Poincaré plot) and calculated its SD1, SD2 and SD1/SD2 defined as the standard deviation of the instantaneous beat-to-beat RRI variability (minor axis of the fitted ellipse), the standard deviation of the continuous long-term RRI variability (major axis of the fitted ellipse) and the axis ratio, respectively¹⁶. They were respectively used as non-linear indexes of (1) the rapid changes in the RRI and thus of the parasympathetic sinus node control¹⁷, (2) the effects of both parasympathetic and sympathetic components on the sinus node activity¹⁸, and (3) the relationship between these components, which is the ratio of the short interval variation to the long interval variation. We also used other non-linear tools to characterize the RRI dynamics, as previously described¹⁰. Recurrence plot analysis quantification (minimal line length (Lmin), mean line length (Lmean), maximal line length (Lmax), divergence (DIV), recurrence rate (REC), determinism (DET) and Shannon entropy (ShanEn)) was performed with embedding dimension, lag, and threshold distance set to m = 10, τ = 1, and r = √mSD of the RRI time series analyzed, respectively. We computed detrended fluctuation analysis α1 and α2 coefficients (alpha1 and alpha2, respectively) with segment length set to n ∈ (4,16) and n ∈ (16,64), respectively. Approximate Entropy (ApEn) and Sample Entropy (SampEn) estimates of each RRI time series were computed with m (embedding dimension) and r (filtering level) set to 2 and 0.2 SD of the RRI time series analyzed, respectively. Finally, we also estimated correlation dimension (CorDim) of the RRI times series¹⁹. Indeed, first, non-linear deterministic measures of heartbeats better detect the autonomic changes related to mortality than the stochastic indexes; and second, what the time-dependent non-linear metrics show as indicators of cardiac vulnerability to lethal arrhythmias are transient non-stationary shifts of dimension of the heartbeat dynamics to a low value²⁰. CorDim quantified the time self-similarity of a signal, i.e., of the RRI time series. CorDim can be considered as an estimation of the lower threshold of the number of degrees of freedom for the underlying system generating the observed data and is typically used as an index of the overall complexity of the system dynamics estimated from a time series²¹. CorDim was computed as described by Grassberger
and Procaccia²² with r = 15% from the standard deviation of the RRI. Accordingly, the first step is to construct the correlation integral (N,m,r) function. The correlation integral counts the fraction of pairs (Xi, Xj) whose distance is smaller than r and is defined as:

$$C(N, m, r) = \frac{1}{N(N-1)} \sum_{i=1}^{N} \sum_{j=i+1}^{N} H(r - |xi - xj|)$$

where Xi and Xj represent phase-space trajectory points, N represents the total amount of phase-space points, and H represents the Heaviside step function, i.e., H(α) = 0 if α < 0 and H(α) = 1 if α ≥ 0. D2 is subsequently computed as:

$$CorDim = \lim_{r \to 0} \frac{\log(C)}{\log(r)}$$

## Outcomes

The main outcome was the ability to express the SSc subtype as a function of HRV, i.e., to classify the SSc clinical form as diffuse or limited of a patient using heart rate dynamics metrics as markers. Secondary outcomes were the link between HRV metrics and SSc skin extent, and between HRV metrics and lung function.

## Statistical analysis

All statistical computations were performed with Python 3.10 LTS.

### Descriptive statistics

The quantitative variables were first tested for normality distribution with the Shapiro-Wilk normality test. Accordingly, means or medians groups' differences (i.e., dcSSc vs lcSSc) were tested performing independent t-tests, paired t-tests, or Wilcoxon-Mann-Whitney tests when necessary for continuous variables, and Wilcoxon-Mann-Whitney or Kruskal-Wallis tests for ordinal variables. The qualitative categorical variables were assessed through modalities headcounts and groups' differences were tested performing Chi-square or Fisher exact tests. Linear correlations between quantitative variables were estimated computing Pearson's and Spearman's correlation coefficients when necessary and hypothesis of correlation was tested performing a t-test. A p-value < 0.05 was considered as statistically significant, nevertheless all variables having an unadjusted univariate p-value < 0.1 (n = 33) were considered for feature engineering. Furthermore, specifically for HRV analysis, and with a dual objective of simplifying the interpretation, particularly in terms of physiological outcomes, and limiting the expansion of type I error, a two-way repeated measures ANOVA analysis with Bonferroni correction was employed. The first factor considered was the clinical form (dcSSc or lcSSc), while the second factor was the position (decubitus or orthostatism). Non-parametric ANOVA (Kruskal-Wallis) was conducted when necessary. For secondary outcomes, simple and multiple linear regressions were performed to predict quantitative variables while simple and multiple logistic regressions were performed to predict bimodal categorical variables. Ordinary Least Squares optimizing method were used in both cases. Significance was set to p-value < 0.05.

### Features engineering

To limit effects of collinearity on features extraction procedure, we identified all pairs of variables being highly correlated (R² > 0.9) and, according to their physiological explainability, excluded the less appropriate variables (n = 9). Quantitative data were then normalized (min-max scaling) for subsequent steps. On one hand, five features selection algorithms (Random Forest, Extra Tree Classifier, Select K Best, Generic Univariate Selection, and Recursive Feature Elimination) were used to rank the 24 selected HRV variables. Ranking was performed on variable importance (explained variance) through a homemade redundancy weighted score taking into account, for each variable, the rank attributed by each algorithm and the dispersion of the rank attributed through the five algorithms. On the other hand, the gamma analysis technique was used to rank the most inter-groups discriminative variables²³,²⁴.

### Models and classification

We used the Synthetic Minority Over-sampling Technique (SMOTE) to manage unbalanced groups' size to avoid any bias in models' learning and classification processes. To limit the number of candidate models and their complexity in order to be compatible with a physiological explainability and a bedside application, we generated only combinations of 1 to 5 variables among the 24 selected HRV variables, i. e. we generated $\sum_{i=1}^{5} C_i^5 = 31$ combinations among the $\sum_{i=1}^{5} C_i^{24} = 55454$ possible combinations. From the 31 combinations of the most discriminative variables generated, 341 classifying models were generated applying on each combination 11 machine learning algorithms that were: logistic regression, decision tree, random forest, linear support vector machine (SVM), radial based function SVM, neuronal network, gradient boosting machine, stacked ensembles, extremely randomized trees, k nearest neighbor, and extreme gradient boosting. Database was randomly split into training (80% of sample size) and validation (20% of sample size) sets. Performances of the 341 classifying models for SSc subtypes classification purposes, i.e., dcSSc or lcSSc, were assessed according to 8 usual metrics that are: mean squared error, logloss, sensitivity, specificity, precision, accuracy, area under curve of the ROC-curves and F1-score. The cross validation was performed with k-fold = 5 and repeated 10 times, providing information on how well a classifier generalizes for test data and further classifications, in particular the expected error range of the classifier thanks to the cross-validation score.
## Ethical approval

The study was approved by the Ethical Committee of Marseille University Hospital and was conducted in accordance with the Declaration of Helsinki. Each subject provided written consent for data analysis. Agreement number is PADS19-352/WZTVF7.

## Results

### Population description

Between January 01, 2018, and December 31, 2021, 60 patients having SSc and assessed for dysautonomia were enrolled. 16 had dcSSc and 42 lcSSc clinical forms according to the internist experts. Two were excluded because of missing data. General characteristics of the studied population and each subgroup are summarized in Table 1. By definition and as expected, the two groups differed clinically mainly in the higher proportion of patients with severe symptoms in the dcSSc group compared with the lcSSc group, as observed through the higher Rodnan (t-statistic = 8.885, df = 20, p < 0.001) and Medsger (U-statistic = 626.5, p < 0.001) scores, and through the presence of sphincter atonia ($\chi^2$-statistic = 26.72, df = 2, p = 0.001) and nonspecific interstitial pneumonia ($\chi^2$-statistic = 44.86, df = 2, p = 0.031). Similarly, esophageal peristaltism dysfunction and diffuse parenchymal lung disease tended to be more frequent in the dcSSc group. As also expected, the two groups biologically differed mainly in the higher proportion of patients with anti-Scl-70 ($\chi^2$-statistic = 209.56, df = 2, p < 0.001) and anti-SSA ($\chi^2$-statistic = 20.5, df = 2, p = 0.001) antibodies in the dcSSc group compared with the lcSSc group that had higher proportion of patients with anti-centromere antibodies ($\chi^2$ = 14.05, df = 2, p = 0.001).

### Lung function

Respiratory parameters measured during lung function test and their interpretation as functional syndromes are summarized in Table 2. As expected, in terms of lung function, the two groups differed only in the higher proportion of patients with restrictive lung disease in the dcSSc group compared with the lcSSc group ($\chi^2$-statistic = 26.74, df = 2, p = 0.001).

### Autonomic influences on sinus node activity

HRV indexes in decubitus then in orthostatism conditions as well as the stand-up induced HRV changes are summarized in Table 3. The two groups differed mainly in a higher HR in the dcSSc group compared with the lcSSc group (mean_HR, clinical form main effect, H-statistic = 11.13, p = 0.003), in decubitus and even more so in orthostatism position (mean_HR, interaction of main effects, H-statistic = 29.93, p < 0.001), HR (mean_HR) being inversely proportional to RRI duration (mean_RR). On the contrary, the differences in observed values for all other HRV markers, especially those derived from spectral and nonlinear analysis, appear to be less systematic. Except for ShanEn, which is similarly higher in the dcSSc group compared with the lcSSc group (ShanEn, clinical form main effect, F-statistic = 6.50, df = 1, p = 0.049) in both decubitus and orthostatic positions (no interaction effect), the other HRV markers differ from one group to another depending on the position considered. This differential effect of verticalization according to patients' group is highlighted through the interaction effect, as seen with LF_power, HF_power, or tot_power. Similarly, a higher length of the mean and the longest diagonal line of the recurrence plot was observed in the dcSSc group compared with the lcSSc group according to position (Lmax, interaction of main effects, H-statistic = 28.92, p < 0.001 and Lmean, interaction of main effects, H-statistic = 13.46, p = 0.011), Lmax being inversely linked to divergence (DIV). Finally, the two groups differed also in terms of HRV parameters changes induced by stand-up. Mainly, mean_HR increase was higher in dcSSc compared to lcSSc groups (t-statistic = 3.034, df = 24, p = 0.006) while all others time domain parameters were lowered during orthostatism and even more lowered in the dcSSc group than in the lcSSc group. Particularly, sd_HR was lowered in the dcSSc group while increased in the lcSSc group (t-statistic = - 3.329, df = 18, p = 0.004). Similarly, SD2 and CorDim were lowered in the dcSSc group while increased in the lcSSc group (t-statistic = - 3.586, df = 18, p = 0.002 and t-statistic = - 3.247, df = 20, p = 0.004, respectively).

### Features importance

Collinearity and correlations between all variables that have been significantly different in the dcSSc group compared to the lcSSc group with a univariate approach are summarized in Figs. 1 and 2. We can note that many HRV indexes are positively or negatively correlated, and mostly of interest, some HRV parameters are correlated to clinical and respiratory quantitative parameters. Particularly, first, mean_HR_S2, SD2_delta, sd_HR_delta, SD2_delta and CorDim_delta are correlated with Rodnan and Medsger scores and with most of the respiratory function parameters, and second, most of the HRV parameters changes (delta) are correlated with clinical scores and lung function parameters. A synthetic view of features importance to describe data variance is proposed in Fig. 3. Features are ranked according to two different metrics, the relative importance of the feature (explained data variance) and the gamma-metric (discriminant power). The first four most efficient features identically identified are SD2_delta, sd_HR_delta, CorDim_delta, and SD2_S2. The fifth important feature is either LF_power_S2 either mean_HR_S2 according to the used metric. As fifth feature we chose mean_HR_S2 because of a more physiological and measuring robustness.

### Models' performances

Table 4 summarizes the performances of the classifying models generated including 1 to 5 selected variables. Concerning sensibility and specificity we can note that models range from 0.87 to 1 and from 0.25 to 0.99 respectively. The most efficient variable taken alone and combining the best sensibility and specificity combination is CorDim with 1 and 0.99 values respectively. But on the other hand, precision, accuracy, and F1-score are the
<table><caption>**Table 1.** General characteristics of the studied population. Values are expressed as mean ± standard deviation or modality size (percentage of the whole population) for quantitative and categorical variables, respectively. The whole studied population (All SSc) are described in the first column and p-values of the comparison between diffuse (dcSSc) and located (IcSSc) cutaneous systemic sclerosis are reported. P < 0.05 was considered as significant and p < 0.1 (bold) was used as threshold for further analysis. BMI: Body mass index, SSc: Systemic sclerosis, NSIP: Nonspecific interstitial pneumonia, IPF: Idiopathic pulmonary fibrosis.</caption><thead><tr><th></th><th>All SSc</th><th>dcSSc</th><th>IcSSc</th><th>p-value</th></tr></thead><tbody><tr><td><strong>Anthropometry</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Age, years</td><td>50.5 ± 13.8</td><td>47.4 ± 11.5</td><td>51.6 ± 14.5</td><td>0.108</td></tr><tr><td>Sex, F/M (%)</td><td>44/14 (76/24)</td><td>9/7 (16/12)</td><td>35/7 (60/12)</td><td><strong>0.068</strong></td></tr><tr><td>Weight, kg</td><td>65.8 ± 11.7</td><td>69.1 ± 12.3</td><td>64.0 ± 10.8</td><td>0.135</td></tr><tr><td>Height, m</td><td>1.66 ± 0.08</td><td>1.71 ± 0.10</td><td>1.65 ± 0.06</td><td><strong>0.044</strong></td></tr><tr><td>BMI, kg/m²</td><td>23.7 ± 3.4</td><td>23.6 ± 3.4</td><td>23.5 ± 3.2</td><td>0.905</td></tr><tr><td><strong>Clinical features and severity</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Raynaud duration, years</td><td>12.9 ± 8.1</td><td>9.5 ± 4.1</td><td>14.2 ± 8.9</td><td><strong>0.01</strong></td></tr><tr><td>SSc duration, years</td><td>10.1 ± 6.3</td><td>8.6 ± 2.8</td><td>10.7 ± 7.1</td><td><strong>0.1</strong></td></tr><tr><td>Rodnan score, n.u.</td><td>6.8 ± 6.4</td><td>15.2 ± 4.8</td><td>3.6 ± 3.2</td><td><strong>&lt;0.001</strong></td></tr><tr><td>Medsger score, n.u.</td><td>4.9 ± 2.5</td><td>7.7 ± 1.6</td><td>3.7 ± 1.8</td><td><strong>&lt;0.001</strong></td></tr><tr><td>Pulmonary hypertension, y/n (%)</td><td>5/53 (9/91)</td><td>2/14 (4/24)</td><td>3/39 (5/67)</td><td>1</td></tr><tr><td>Oesophagal aperistaltism, y/n (%)</td><td>31/27 (53/47)</td><td>13/3 (22/5)</td><td>18/24 (31/42)</td><td><strong>0.064</strong></td></tr><tr><td>Sphincter atonia, y/n (%)</td><td>22/36 (28/72)</td><td>13/3 (22/5)</td><td>9/33 (16/57)</td><td><strong>&lt;0.001</strong></td></tr><tr><td><strong>Lung disease</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Dyspnea, I/II/III/IV NYHA scale</td><td>33/23/2/0</td><td>7/7/2/0</td><td>26/16/0/0</td><td>0.176</td></tr><tr><td>Thorax TDM, y/n (%)</td><td>43/15 (74/26)</td><td>12/4 (21/7)</td><td>31/11 (53/19)</td><td>0.984</td></tr><tr><td>Diffuse parenchymal lung disease, y/n (%)</td><td>13/30 (30/70)</td><td>6/6 (14/14)</td><td>7/24 (16/56)</td><td><strong>0.076</strong></td></tr><tr><td>NSIP, y/n (%)</td><td>3/40 (7/93)</td><td>2/10 (5/23)</td><td>1/30 (2/70)</td><td><strong>0.031</strong></td></tr><tr><td>IPF, y/n (%)</td><td>3/40 (7/93)</td><td>1/11 (2/25)</td><td>2/29 (5/68)</td><td>0.965</td></tr><tr><td>Fibrosis, y/n (%)</td><td>10/33 (23/77)</td><td>4/8 (9/19)</td><td>6/25 (14/58)</td><td>0.472</td></tr><tr><td>Fibrosis &gt; 20%, y/n (%)</td><td>8/35 (19/81)</td><td>3/9 (7/21)</td><td>5/26 (12/60)</td><td>0.705</td></tr><tr><td><strong>Biology</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Anti-Scl-70, y/n (%)</td><td>31/26 (54/46)</td><td>15/1 (26/2)</td><td>16/25 (28/44)</td><td><strong>&lt;0.001</strong></td></tr><tr><td>Anti-centromere, y/n (%)</td><td>24/33 (42/58)</td><td>0/16 (0/28)</td><td>24/17 (42/30)</td><td><strong>0.001</strong></td></tr><tr><td>Anti-SSA, y/n (%)</td><td>8/49 (14/86)</td><td>8/8 (14/14)</td><td>0/41 (0/72)</td><td><strong>0.001</strong></td></tr><tr><td>Anti-SSB, y/n (%)</td><td>5/52 (9/91)</td><td>5/11 (9/19)</td><td>0/41 (0/72)</td><td><strong>0.054</strong></td></tr><tr><td>Anti-U1 RNP, y/n (%)</td><td>0/57 (0/100)</td><td>0/16 (0/28)</td><td>0/41 (0/72)</td><td>1</td></tr><tr><td>Anti-RNA polymerase III, y/n (%)</td><td>1/56 (2/98)</td><td>1/15 (2/26)</td><td>0/41 (0/72)</td><td>0.918</td></tr><tr><td>Anti-PM/Scl y/n (%)</td><td>0/57 (0/100)</td><td>0/16 (0/28)</td><td>0/41 (0/72)</td><td>1</td></tr><tr><td>Anti-fibrillarin, y/n (%)</td><td>1/56 (2/98)</td><td>0/16 (0/28)</td><td>1/40 (2/70)</td><td>0.988</td></tr><tr><td><strong>Treatments</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Calcium channel blockers, y/n (%)</td><td>35/23 (60/40)</td><td>7/9 (12/16)</td><td>28/14 (48/24)</td><td><strong>0.011</strong></td></tr><tr><td>Diltiazem, y/n (%)</td><td>5/53 (9/91)</td><td>2/14 (3/24)</td><td>3/39 (5/68)</td><td>0.576</td></tr><tr><td>Non-diltiazem, y/n (%)</td><td>30/28 (52/48)</td><td>5/11 (9/19)</td><td>25/17 (43/29)</td><td><strong>&lt;0.001</strong></td></tr><tr><td>Proton-pump inhibitors, y/n (%)</td><td>10/48 (17/83)</td><td>4/12 (7/21)</td><td>6/36 (10/62)</td><td>0.276</td></tr><tr><td>Steroids, y/n (%)</td><td>2/56 (3/97)</td><td>1/15 (2/26)</td><td>1/42 (2/70)</td><td>0.585</td></tr><tr><td>L-thyroxin, y/n (%)</td><td>6/52 (10/90)</td><td>0/16 (0/28)</td><td>6/36 (10/62)</td><td>0.651</td></tr><tr><td>Chloroquin, y/n (%)</td><td>4/54 (7/93)</td><td>0/16 (0/28)</td><td>4/38 (7/65)</td><td>0.827</td></tr><tr><td>Sildenafil, y/n (%)</td><td>7/51 (12/88)</td><td>3/13 (5/22)</td><td>4/38 (7/66)</td><td>0.309</td></tr><tr><td>Mycophenolate, y/n (%)</td><td>12/46 (21/79)</td><td>9/7 (16/12)</td><td>3/39 (5/67)</td><td><strong>&lt;0.001</strong></td></tr><tr><td>Bosentan, y/n (%)</td><td>8/50 (14/86)</td><td>3/13 (5/22)</td><td>5/37 (9/64)</td><td>0.524</td></tr><tr><td>Other drugs, y/n (%)</td><td>14/44 (24/76)</td><td>3/13 (5/22)</td><td>11/31 (19/54)</td><td>0.466</td></tr></tbody></table>

lowest of all tested models (0.749, 0.771, and 0.823 respectively). Besides, it should be noted that the proposed models, except for the one including only CorDim, appear to have a rather average sensitivity ranging from 0.25 to 0.81. Moreover and interestingly, we can note that F1-score that is an integrative metrics (harmonic mean of precision and accuracy) is not a linear function of the number of included variables. F1-score ranges from 0.823 for the best 1-variable model to a maximum of 0.947 for the 4-variables best model while 0.946 is reached by the 3-variables best model. If specificity and/or precision is an objective we can note that best models are two 3-variables and two 4-variables models including SD2_delta, SD2_S2, sd_HR_delta, and CorDim_delta.
<table><caption>**Table 2.** Lung function characteristics of the studied population. Values are expressed as mean ± standard deviation or modality size (percentage of the whole population) for quantitative and categorical variables, respectively. The whole studied population (All SSc) are described in the first column and p-values of the comparison between diffuse (dcSSc) and located (lcSSc) cutaneous systemic sclerosis are reported. P < 0.05 was considered as significant and p < 0.1 (bold) was used as threshold for further analysis. TLC, total lung capacity; FRC, functional residual capacity; slow VC, slow vital capacity; forced VC, forced vital capacity; RV/ TLC, residual volume/total lung capacity ratio; FEV1, first 1-second forced expiratory volume; FEV1/forced VC, modified Tiffeneau-Pinelli index; MMEF, Maximal mid-expiratory flow; Raw, central airway resistance; DLCOsb, single-breath carbon monoxide diffusion capacity; KCOc, corrected carbon monoxide transfer coefficient; SpO2, oxyhemoglobin saturation rate estimated from pulse oximetry.</caption><thead><tr><th></th><th>All SSc</th><th>dcSSc</th><th>lcSSc</th><th>p-value</th></tr></thead><tbody><tr><td colspan="5"><strong>Lung function, values (n=58)</strong></td></tr><tr><td>TLC, L</td><td>4.87 ± 1.11</td><td>4.10 ± 0.89</td><td>5.15 ± 1.06</td><td><strong>&lt;0.001</strong></td></tr><tr><td>% of predicted</td><td>91 ± 21</td><td>71 ± 17</td><td>98 ± 17</td><td><strong>&lt;0.001</strong></td></tr><tr><td>FRC, L (%)</td><td>2.83 ± 0.60</td><td>2.53 ± 0.44</td><td>2.95 ± 0.62</td><td><strong>0.007</strong></td></tr><tr><td>% of predicted</td><td>99 ± 21</td><td>83 ± 14</td><td>105 ± 21</td><td><strong>&lt;0.001</strong></td></tr><tr><td>Slow VC, L (%)</td><td>3.03 ± 0.85</td><td>2.43 ± 0.63</td><td>3.26 ± 0.81</td><td><strong>&lt;0.001</strong></td></tr><tr><td>% of predicted</td><td>91 ± 27</td><td>65 ± 24</td><td>101 ± 20</td><td><strong>&lt;0.001</strong></td></tr><tr><td>Forced VC, L (%)</td><td>2.94 ± 0.87</td><td>2.31 ± 0.63</td><td>3.18 ± 0.83</td><td><strong>&lt;0.001</strong></td></tr><tr><td>% of predicted</td><td>91 ± 28</td><td>64 ± 26</td><td>101 ± 20</td><td><strong>&lt;0.001</strong></td></tr><tr><td>RV/TLC, no unit (%)</td><td>37.75 ± 7.88</td><td>40.63 ± 7.64</td><td>36.72 ± 7.79</td><td><strong>0.084</strong></td></tr><tr><td>% of predicted</td><td>105 ± 22</td><td>120 ± 26</td><td>100 ± 18</td><td><strong>&lt;0.001</strong></td></tr><tr><td>FEV<sub>1</sub>, L (%)</td><td>2.37 ± 0.67</td><td>1.90 ± 0.42</td><td>2.54 ± 0.67</td><td><strong>&lt;0.001</strong></td></tr><tr><td>% of predicted</td><td>87 ± 25</td><td>63 ± 21</td><td>96 ± 19</td><td><strong>&lt;0.001</strong></td></tr><tr><td>FEV<sub>1</sub>/forced VC, %, (%)</td><td>81.41 ± 7.70</td><td>83.75 ± 8.88</td><td>80.51 ± 7.12</td><td>0.163</td></tr><tr><td>% of predicted</td><td>103 ± 9</td><td>105 ± 10</td><td>102 ± 8</td><td>0.283</td></tr><tr><td>MMEF, L (%)</td><td>2.47 ± 1.06</td><td>2.12 ± 0.65</td><td>2.6 ± 1.17</td><td>0.443</td></tr><tr><td>% of predicted</td><td>73 ± 28</td><td>58 ± 16</td><td>79 ± 29</td><td><strong>0.001</strong></td></tr><tr><td>Raw cmH<sub>2</sub>O*s/L (%)</td><td>1.90 ± 0.74</td><td>1.91 ± 0.63</td><td>1.90 ± 0.78</td><td>0.711</td></tr><tr><td>% of predicted</td><td>62 ± 24</td><td>63 ± 21</td><td>62 ± 26</td><td>0.711</td></tr><tr><td>DLCO<sub>sb</sub>, mL/mmHg/min (%)</td><td>16.38 ± 5.19</td><td>14.48 ± 3.14</td><td>17.08 ± 5.63</td><td><strong>0.036</strong></td></tr><tr><td>% of predicted</td><td>65 ± 18</td><td>53 ± 12</td><td>69 ± 18</td><td><strong>&lt;0.001</strong></td></tr><tr><td>KCO<sub>c</sub>, mL/mmHg/min/L (%)</td><td>3.76 ± 0.78</td><td>3.99 ± 0.78</td><td>3.69 ± 0.77</td><td>0.225</td></tr><tr><td>% of predicted</td><td>80 ± 16</td><td>86 ± 16</td><td>78 ± 15</td><td>0.150</td></tr><tr><td>SpO<sub>2</sub>, %</td><td>99 ± 0.99</td><td>98 ± 1</td><td>99 ± 1</td><td>0.5</td></tr><tr><td colspan="5"><strong>Lung function, interpretation (n=58)</strong></td></tr><tr><td>Restrictive lung disease, y/n (%)</td><td>16/42 (28/72)</td><td>9/7 (16/12)</td><td>7/35 (12/60)</td><td><strong>&lt;0.001</strong></td></tr><tr><td>Lung distension, y/n (%)</td><td>1/57 (2/98)</td><td>0/16 (0/28)</td><td>1/41 (2/70)</td><td>0.988</td></tr><tr><td>Obstructive lung, y/n (%)</td><td>5/53 (9/91)</td><td>1/15 (2/26)</td><td>4/38 (7/65)</td><td>0.681</td></tr><tr><td>Lung diffusion impairment, y/n (%)</td><td>11/47 (19/81)</td><td>2/14 (3/24)</td><td>9/33 (16/57)</td><td>0.216</td></tr></tbody></table>

## Secondary outcomes

The link between HRV metrics and SSc skin extent appears through some significant regressions, the most significant simple linear regressions being Rodnan score = -165.063 × SD2_delta + 6.261 (R² = 0.268, F-statistic: 19.80, df = 1, p < 0.001) and Medsger score = -0.062 × LF_power_prc_delta + 5.143 (R² = 0.210, F-statistic: 14.36, df = 1, p < 0.001). The most significant multiple linear regressions provide Rodnan score estimation from mean_RR_S2, mean_HR_S2, ApEn_S2, alpha1_S2, Lmean_S2, DIV_S2, DET_S2, and ShanEn_S2 (R² = 0.842, F-statistic: 4.453, df = 30, p < 0.001) and Medsger score estimation from alpha1_S2, CorDim_S2, and Lmean_S2 (R² = 0.756, F-statistic: 2.393, df = 30, p = 0.015). Similarly, the link between HRV metrics and lung function appears through some significant regressions, the most significant simple logistic regression expressing the presence of a restrictive syndrome as a function of mean_HR_S2 (pseudo R² = 0.277, likelihood statistic: -32.542, df = 1, p < 0.001) and the most significant multiple logistic regression expressing the presence of a restrictive syndrome as a function of SD2_S2 and sd_RR_S2 (pseudo R² = 0.55, likelihood statistic: -32.542, df = 8, p < 0.001).

## Discussion

### Major findings

We aimed to provide classifying models of SSc subtypes, i.e., models that distinguish dcSSc from lcSSc clinical forms, from HRV linear and nonlinear metrics, and easy to compute at the patient's bedside. As expected, at rest in decubitus condition, only HR was clearly different in SSc subgroups, HR being higher in dcSSc. Transition
<table><thead><tr><th rowspan="2"></th><th colspan="2">Decubitus</th><th colspan="2">Othostatism</th><th>Clinical form</th><th>Position</th><th>Interactions</th><th colspan="3">Stand-up-induced changes</th></tr><tr><th>dcSSc</th><th>lcSSc</th><th>dcSSc</th><th>lcSSc</th><th>p value</th><th>p value</th><th>p value</th><th>DcSSc</th><th>lcSSc</th><th>p value</th></tr></thead><tbody><tr><td colspan="11"><strong>Time domain</strong></td></tr><tr><td>mean_RR (ms)</td><td>825 ± 120</td><td>894 ± 113</td><td>698 ± 117</td><td>795 ± 102</td><td><strong>0.002*</strong></td><td>&lt;0.001</td><td>1</td><td>-127 ± 48</td><td>-99 ± 52</td><td>0.072*</td></tr><tr><td>sd_RR (ms)</td><td>40 ± 22</td><td>33 ± 14</td><td>22 ± 11</td><td>34 ± 14</td><td>0.590</td><td>0.555</td><td><strong>0.035*</strong></td><td>-18 ± 22</td><td>1 ± 9</td><td><strong>0.005*</strong></td></tr><tr><td>mean_HR (bpm)</td><td>74.0 ± 10.0</td><td>68.0 ± 9.0</td><td>88.0 ± 13.0</td><td>77.0 ± 11.0</td><td><strong>0.003*</strong></td><td>&lt;0.001</td><td><strong>&lt;0.001</strong></td><td>14.0 ± 6.0</td><td>9.0 ± 5.0</td><td><strong>0.006*</strong></td></tr><tr><td>sd_HR (bpm)</td><td>4.0 ± 2.0</td><td>3.0 ± 1.0</td><td>3.0 ± 1.0</td><td>3.0 ± 1.0</td><td>1</td><td>0.426</td><td><strong>0.047*</strong></td><td>-1.0 ± 2.0</td><td>1.0 ± 1.0</td><td><strong>0.004*</strong></td></tr><tr><td>RMSSD (ms)</td><td>35 ± 29</td><td>28 ± 18</td><td>14 ± 13</td><td>22 ± 16</td><td>0.890</td><td>&lt;0.001</td><td><strong>0.002*</strong></td><td>-21 ± 29</td><td>-6 ± 10</td><td>0.052*</td></tr><tr><td>NN50 (count)</td><td>53.0 ± 61.0</td><td>30.0 ± 52.0</td><td>12.0 ± 30.0</td><td>22.0 ± 44.0</td><td>1</td><td>0.046</td><td>0.068</td><td>-41.0 ± 60.0</td><td>-8.0 ± 26.0</td><td><strong>0.048*</strong></td></tr><tr><td>pNN50 (%)</td><td>15.0 ± 17.0</td><td>9.0 ± 16.0</td><td>4.0 ± 9.0</td><td>6.0 ± 12.0</td><td>1</td><td>0.030</td><td><strong>0.048</strong></td><td>-12.0 ± 16.0</td><td>-3.0 ± 8.0</td><td>0.054*</td></tr><tr><td colspan="11"><strong>Frequency domain</strong></td></tr><tr><td>LF_peak (Hz)</td><td>0.074 ± 0.022</td><td>0.076 ± 0.028</td><td>0.074 ± 0.025</td><td>0.073 ± 0.027</td><td>1</td><td>1</td><td>1</td><td>-0.000 ± 0.00</td><td>-0.003 ± 0.03</td><td>0.648</td></tr><tr><td>HF_peak (Hz)</td><td>0.264 ± 0.08</td><td>0.257 ± 0.07</td><td>0.24 ± 0.071</td><td>0.214 ± 0.069</td><td>1</td><td>0.008</td><td><strong>0.038</strong></td><td>-0.024 ± 0.101</td><td>-0.044 ± 0.105</td><td>0.393</td></tr><tr><td>LF_power (ms²)</td><td>750 ± 1574</td><td>369 ± 508</td><td>155 ± 218</td><td>434 ± 434</td><td>0.142</td><td>1</td><td><strong>0.013*</strong></td><td>-595 ± 1602</td><td>65 ± 464</td><td>0.124</td></tr><tr><td>LF_power_prc (%)</td><td>31.0 ± 12.0</td><td>30.0 ± 13.0</td><td>27.0 ± 17.0</td><td>37.0 ± 16.0</td><td>0.461</td><td>1</td><td>0.302*</td><td>-4.0 ± 16.0</td><td>7.0 ± 19.0</td><td><strong>0.040*</strong></td></tr><tr><td>LF_power_nu (n.u.)</td><td>54.0 ± 18.0</td><td>53.0 ± 18.0</td><td>71.0 ± 24.0</td><td>71.0 ± 19.0</td><td>1</td><td>&lt;0.001</td><td><strong>&lt;0.001</strong></td><td>17.0 ± 20.0</td><td>18.0 ± 23.0</td><td>0.485</td></tr><tr><td>HF_power (ms²)</td><td>712 ± 1385</td><td>326 ± 377</td><td>89 ± 152</td><td>192 ± 265</td><td>0.946</td><td>&lt;0.001</td><td><strong>&lt;0.001*</strong></td><td>-623 ± 1391</td><td>-134 ± 289</td><td>0.183</td></tr><tr><td>HF_power_prc (%)</td><td>27.0 ± 13.0</td><td>27.0 ± 17.0</td><td>13.0 ± 16.0</td><td>16.0 ± 15.0</td><td>1</td><td>&lt;0.001</td><td><strong>&lt;0.001</strong></td><td>-13.0 ± 13.0</td><td>-11.0 ± 16.0</td><td>0.686</td></tr><tr><td>HF_power_nu (n.u.)</td><td>46.0 ± 18.0</td><td>46.0 ± 18.0</td><td>29.0 ± 24.0</td><td>28.0 ± 19.0</td><td>1</td><td>&lt;0.001</td><td><strong>&lt;0.001</strong></td><td>-17.0 ± 20.0</td><td>-18.0 ± 23.0</td><td>0.496</td></tr><tr><td>LF_HF_power</td><td>1.54 ± 1.08</td><td>1.61 ± 1.55</td><td>4.72 ± 4.05</td><td>4.59 ± 4.6</td><td>1</td><td>&lt;0.001</td><td><strong>&lt;0.001</strong></td><td>3.18 ± 3.67</td><td>2.98 ± 4.9</td><td>0.866</td></tr><tr><td>tot_power (ms²)</td><td>2055 ± 3373</td><td>1135 ± 1004</td><td>559 ± 567</td><td>1218 ± 1081</td><td>0.657</td><td>0.502</td><td>0.092*</td><td>-1496 ± 3474</td><td>83 ± 982</td><td>0.092*</td></tr><tr><td colspan="11"><strong>Nonlinear domain</strong></td></tr><tr><td>SD1 (ms)</td><td>25 ± 20</td><td>20 ± 13</td><td>10 ± 9</td><td>15 ± 11</td><td>0.890</td><td>&lt;0.001</td><td><strong>0.002*</strong></td><td>-15 ± 20</td><td>-4 ± 7</td><td>0.052*</td></tr><tr><td>SD2 (ms)</td><td>50 ± 24</td><td>42 ± 17</td><td>29 ± 13</td><td>45 ± 19</td><td>0.898</td><td>0.559</td><td><strong>0.034*</strong></td><td>-21 ± 26</td><td>3 ± 13</td><td><strong>0.002*</strong></td></tr><tr><td>ApEn</td><td>1.139 ± 0.067</td><td>1.112 ± 0.11</td><td>1.006 ± 0.176</td><td>1.071 ± 0.105</td><td>1</td><td>0.003</td><td><strong>0.015</strong></td><td>-0.133 ± 0.169</td><td>-0.041 ± 0.138</td><td><strong>0.034*</strong></td></tr><tr><td>SampEn</td><td>1.561 ± 0.275</td><td>1.58 ± 0.268</td><td>1.171 ± 0.359</td><td>1.298 ± 0.25</td><td>0.832</td><td>&lt;0.001</td><td>1</td><td>-0.391 ± 0.415</td><td>-0.281 ± 0.282</td><td>0.257</td></tr><tr><td>alpha1</td><td>1.045 ± 0.288</td><td>1.047 ± 0.248</td><td>1.261 ± 0.355</td><td>1.308 ± 0.313</td><td>1</td><td>&lt;0.001</td><td><strong>&lt;0.001</strong></td><td>0.216 ± 0.319</td><td>0.261 ± 0.285</td><td>0.636</td></tr><tr><td>alpha2</td><td>0.964 ± 0.161</td><td>0.928 ± 0.177</td><td>1.12 ± 0.206</td><td>0.979 ± 0.221</td><td>0.133</td><td>0.12</td><td>0.801*</td><td>0.156 ± 0.215</td><td>0.051 ± 0.239</td><td>0.054*</td></tr><tr><td>CorDim</td><td>2.082 ± 1.703</td><td>1.358 ± 1.367</td><td>0.759 ± 1.163</td><td>1.449 ± 1.367</td><td>1</td><td>0.745</td><td>0.115*</td><td>-1.323 ± 1.619</td><td>0.091 ± 1.039</td><td><strong>0.004*</strong></td></tr><tr><td>Lmin (beats)</td><td>2.0 ± 0.0</td><td>2.0 ± 0.0</td><td>2.0 ± 0.0</td><td>2.0 ± 0.0</td><td>1</td><td>1</td><td>1</td><td>0.0 ± 0.0</td><td>0.0 ± 0.0</td><td>1.0</td></tr><tr><td>Lmax (beats)</td><td>235.0 ± 126.0</td><td>159.0 ± 75.0</td><td>320.0 ± 145.0</td><td>285.0 ± 115.0</td><td>0.194</td><td>&lt;0.001</td><td><strong>&lt;0.001*</strong></td><td>85.0 ± 115.0</td><td>126.0 ± 108.0</td><td>0.168</td></tr><tr><td>Lmean (beats)</td><td>12.5 ± 3.7</td><td>11.0 ± 4.0</td><td>15.0 ± 4.7</td><td>12.6 ± 3.5</td><td>0.090</td><td>0.009</td><td><strong>0.011*</strong></td><td>2.5 ± 6.0</td><td>1.6 ± 3.5</td><td>0.496</td></tr><tr><td>DIV</td><td>0.006 ± 0.003</td><td>0.008 ± 0.006</td><td>0.004 ± 0.003</td><td>0.005 ± 0.005</td><td>0.194</td><td>&lt;0.001</td><td><strong>&lt;0.001*</strong></td><td>-0.002 ± 0.003</td><td>-0.003 ± 0.005</td><td>0.163</td></tr><tr><td>REC (%)</td><td>0.368 ± 0.118</td><td>0.315 ± 0.097</td><td>0.399 ± 0.097</td><td>0.373 ± 0.084</td><td>0.198</td><td>0.021</td><td>1</td><td>0.032 ± 0.139</td><td>0.058 ± 0.092</td><td>0.686</td></tr><tr><td>DET (%)</td><td>98.0 ± 1.5</td><td>97.6 ± 1.3</td><td>98.9 ± 1</td><td>98.7 ± 1.2</td><td>0.200</td><td>&lt;0.001</td><td><strong>&lt;0.001</strong></td><td>0.9 ± 1.3</td><td>1.1 ± 1.1</td><td>0.374</td></tr><tr><td>ShanEn</td><td>3.262 ± 0.312</td><td>3.111 ± 0.302</td><td>3.475 ± 0.357</td><td>3.293 ± 0.31</td><td><strong>0.049*</strong></td><td>0.006</td><td>1</td><td>0.212 ± 0.44</td><td>0.183 ± 0.265</td><td>0.776</td></tr></tbody></table>

**Table 3.** Heart rate variability of the studied population. Values are expressed as mean ± standard deviation. Diffuse (dcSSc) and located (lcSSc) cutaneous systemic sclerosis are described for decubitus and orthostatism and adjusted (Bonferroni's correction) p-values of the clinical form and position main effects as their interaction are reported. Stand-up-induced changes of each group are also reported. P < 0.05 was considered as significant (bold). Unadjusted p-values < 0.1 from univariate analysis (\* mark) was used as threshold for further analysis. Time domain variables are: mean_RR, mean RR intervals duration; sd_RR, standard deviation of RR intervals duration; mean_HR, mean heart rate; sd_HR, standard deviation of heart rate; RMSSD, root mean square of successive differences of RR intervals duration; NN50, count of successive RR intervals duration differing from more than 50 ms; pNN50, percentage of RR intervals duration differing from more than 50 ms. Frequency domain variables are: LF_peak and HF_peak, low frequency (0.04–0.15 Hz) and high frequency (0.15–0.40 Hz) having the highest power peaks, respectively; LF_power, low frequency power; LF_power_prc, percentage of low frequency power; LF_power_nu, normalized low frequency power; HF_power, high frequency power; HF_power_prc, percentage of high frequency power; HF_power_nu, normalized high frequency power; LF_HF_power, low frequency/high frequency power ratio; tot_power, total spectral power. Nonlinear variables are: SD1, minor axis of the fitted ellipse on RR intervals Poincaré plot; SD2, major axis of the fitted ellipse on RR intervals Poincaré plot; ApEn, approximate entropy; SamEn, sample entropy; alpha1, first (short-term) alpha coefficient of detrended fluctuation analysis of RR intervals time series; alpha2, second (long-term) alpha coefficient of detrended fluctuation analysis of RR intervals time series; CorDim, correlation dimension; Lmin, minimal line length; Lmax, maximal line length; Lmean, mean line length; DIC, divergence; REC, recurrence rate; DET, determinism; ShanEn, Shannon entropy.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Correlations between heart rate variability metrics. The HRV metrics that were significantly different (*p* < 0.05) in the dcSSc and lcSSc subgroups were tested for correlation and are represented in the correlation matrix. Correlations between HRV metrics from time-, frequency-, and nonlinear- domains measured in supine (_S1_) as orthostatism (_S2_) conditions and HRV metrics stand-up induced changes (_delta_) are reported. R² is color-coded from -1 (dark red) to +1 (dark blue) when correlation was significant (*p* < 0.05).</figcaption></figure>

from decubitus to orthostatism increased HR in all SSc patients. HR reached higher magnitude in dcSSc and was associated with a lower overall variability (tot_power). Stand-up-induced HRV changes revealed differential effects of orthostatism in SSc subtypes, effects that were mainly higher HR increase associated with (1) loss of time-domain variability and (2) loss of nonlinear proprieties of HRV (sd_RR, sd_HR, SD1 and SD2, ApEn, CorDim, and recurrence plot quantifying indexes) in dcSSc compared to lcSSc patients. Accordingly, efficient classifying models of SSc subtypes have been provided using from 1 to 4 HRV metrics as input variables reaching F1-score of overall performance up to almost 0.95. To resume, it's possible to distinguish dcSSc from lcSSc clinical forms at the patient's bedside with a simple digital ECG analysis.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Correlations between clinical features, heart rate variability and lung function metrics. Correlations between HRV metrics from time-, frequency-, and nonlinear- domains measured in supine (S1) as orthostatism (S2) conditions and HRV metrics stand-up induced changes (delta), and clinical features as well lung function metrics are reported. R² is color-coded from -1 (dark red) to +1 (dark blue) when correlation was significant (p < 0.05).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Features importance. The features importance of heart rate variability metrics is represented according to two different features importance ranking techniques: relative importance (top histograms), and gamma-metric (bottom histograms). Features importance of HRV metrics from time-, frequency-, and nonlinear- domains measured in supine (_S1_) as orthostatism (_S2_) conditions and HRV metrics stand-up induced changes (_delta_) was assessed. The four most important features (red arrows) were: SD2_delta, stand-up induced change of the major axis length of the fitted ellipse on RR intervals Poincaré plot; sd_HR_delta, stand-up induced change of the heart rate standard deviation; CorDim_delta, stand-up induced change of the correlation dimension; SD2_S2, major axis length of the fitted ellipse on RR intervals Poincaré plot in orthostatism condition.</figcaption></figure>

## Autonomic influences on sinus node activity

Time domain and frequency domain HRV metrics have been used to investigate the cardiovascular consequences of SSc since 1994². On one hand, SSc is known to alter autonomic nervous system and its cardiovascular component²⁵ and on the other hand HRV has been considered to inform on autonomic nervous system modulation of sinus node rhythmic activity¹³. High frequency (around 0.25 Hz) being related to respiration via the parasympathetics and low frequency (around 0.1 Hz) to a sympathetic Eigen-oscillation of the baroreflex to blood pressure and heart rate system²⁶. Typically, SSc is described as increasing resting HR through a predominant sympathetic activity and vagal withdrawal compared to healthy controls⁹,²⁷⁻²⁹ and as blunting cardiovascular response to orthostatism⁹,²⁹. The HR and HRV we observed in dcSSc and lcSSc subtypes at rest and during orthostatism does not evidence a so clear phenotype-dependent sympathetic activation neither parasympathetic withdrawal. In decubitus position, when comparing dcSSc and lcSSc groups using multiple univariate unpaired t-tests as well as ANOVA, we found no evidence for any difference in autonomic modulation of sinus activity, as commonly accepted¹³. For 5-min RR intervals time series this usually includes some HRV spectral indexes that are LF_power (ms²), LF_power_prc (%), LF_power_nu (n.u.), HF_power (ms²), HF_power_prc (%), HF_power_nu (n.u.) and LF_HF_power. On the contrary, in dcSSc compared to lcSSc,
<table><caption>**Table 4.** Performances of classifying models. SSc subtypes classifying models including 1 to 5 most informative and discriminant HRV variables are presented. Variables used are: SD2_delta, stand-up induced change of the major axis length of the fitted ellipse on RR intervals Poincaré plot; sd_HR_delta, stand-up induced change of the heart rate standard deviation; CorDim_delta, stand-up induced change of the correlation dimension; SD2_S2, major axis length of the fitted ellipse on RR intervals Poincaré plot in orthostatism condition; mean_HR_S2, mean heart rate in orthostatism condition. Classifier algorithms used were neuronal networks; stacked ensembles; and GBM, gradient boosting machines. Models' performances were assessed through MSE, mean squared error; LogLoss, logarithmic loss; Sensibility; Specificity; Precision; Accuracy; AUC, area under curve; and F1-score.</caption><thead><tr><th>Model variables</th><th>Classifier algorithm</th><th>MSE</th><th>LogLoss</th><th>Sensibility</th><th>Specificity</th><th>Precision</th><th>Accuracy</th><th>AUC</th><th>F1-score</th></tr></thead><tbody><tr><td colspan="10"><strong>1 variable</strong></td></tr><tr><td>mean_HR_S2</td><td>Neuronal network</td><td>0.17</td><td>0.578</td><td>0.96</td><td>0.583</td><td>0.858</td><td>0.869</td><td>0.762</td><td>0.897</td></tr><tr><td>SD2_delta</td><td>Neuronal network</td><td>0.157</td><td>0.433</td><td>0.95</td><td>0.5</td><td>0.841</td><td>0.847</td><td>0.875</td><td>0.88</td></tr><tr><td>SD2_S2</td><td>Neuronal network</td><td>0.162</td><td>0.492</td><td>0.92</td><td>0.563</td><td>0.837</td><td>0.829</td><td>0.782</td><td>0.852</td></tr><tr><td>sd_HR_delta</td><td>Neuronal network</td><td>0.182</td><td>1.192</td><td>1</td><td>0.25</td><td>0.72</td><td>0.762</td><td>0.759</td><td>0.828</td></tr><tr><td>CorDim_delta</td><td>Stacked Ensemble</td><td>0.188</td><td>0.53</td><td>1</td><td>1</td><td>0.729</td><td>0.749</td><td>0.771</td><td>0.823</td></tr><tr><td colspan="10"><strong>2 variables</strong></td></tr><tr><td>SD2_delta, mean_HR_S2</td><td>GBM</td><td>0.117</td><td>0.341</td><td>1</td><td>0.583</td><td>0.858</td><td>0.891</td><td>0.9</td><td>0.919</td></tr><tr><td>CorDim_delta, mean_HR_S2</td><td>Neuronal network</td><td>0.154</td><td>0.465</td><td>1</td><td>0.583</td><td>0.858</td><td>0.891</td><td>0.848</td><td>0.919</td></tr><tr><td>SD2_delta, CorDim_delta</td><td>GBM</td><td>0.155</td><td>0.468</td><td>0.92</td><td>0.708</td><td>0.916</td><td>0.891</td><td>0.88</td><td>0.903</td></tr><tr><td>sd_HR_delta, mean_HR_S2</td><td>GBM</td><td>0.137</td><td>0.423</td><td>1</td><td>0.521</td><td>0.825</td><td>0.869</td><td>0.888</td><td>0.901</td></tr><tr><td>SD2_S2, mean_HR_S2</td><td>Neuronal network</td><td>0.159</td><td>0.453</td><td>1</td><td>0.604</td><td>0.833</td><td>0.871</td><td>0.88</td><td>0.901</td></tr><tr><td>sd_HR_delta, CorDim_delta</td><td>Neuronal network</td><td>0.186</td><td>0.546</td><td>0.91</td><td>0.625</td><td>0.882</td><td>0.869</td><td>0.88</td><td>0.888</td></tr><tr><td>sd_HR_delta, SD2_S2</td><td>Neuronal network</td><td>0.155</td><td>0.503</td><td>0.87</td><td>0.688</td><td>0.922</td><td>0.869</td><td>0.83</td><td>0.878</td></tr><tr><td>SD2_delta, SD2_S2</td><td>Neuronal network</td><td>0.141</td><td>0.438</td><td>1</td><td>0.396</td><td>0.783</td><td>0.824</td><td>0.859</td><td>0.873</td></tr><tr><td>CorDim_delta, SD2_S2</td><td>Neuronal network</td><td>0.199</td><td>0.616</td><td>0.95</td><td>0.375</td><td>0.806</td><td>0.802</td><td>0.802</td><td>0.854</td></tr><tr><td>SD2_delta, sd_HR_delta</td><td>Neuronal network</td><td>0.162</td><td>0.515</td><td>1</td><td>0.313</td><td>0.738</td><td>0.784</td><td>0.811</td><td>0.841</td></tr><tr><td colspan="10"><strong>3 variables</strong></td></tr><tr><td>SD2_delta, SD2_S2, mean_HR_S2</td><td>GBM</td><td>0.101</td><td>0.311</td><td>1</td><td>0.771</td><td>0.902</td><td>0.936</td><td>0.943</td><td>0.946</td></tr><tr><td>SD2_delta, CorDim_delta, mean_HR_S2</td><td>Neuronal network</td><td>0.126</td><td>0.434</td><td>0.96</td><td>0.688</td><td>0.916</td><td>0.911</td><td>0.903</td><td>0.935</td></tr><tr><td>CorDim_delta, SD2_S2, mean_HR_S2</td><td>Neuronal network</td><td>0.133</td><td>0.54</td><td>1</td><td>0.563</td><td>0.881</td><td>0.889</td><td>0.826</td><td>0.929</td></tr><tr><td>sd_HR_delta, CorDim_delta, mean_HR_S2</td><td>Neuronal network</td><td>0.123</td><td>0.465</td><td>0.95</td><td>0.625</td><td>0.898</td><td>0.891</td><td>0.868</td><td>0.913</td></tr><tr><td>SD2_delta, sd_HR_delta, SD2_S2</td><td>Neuronal network</td><td>0.165</td><td>0.515</td><td>1</td><td>0.5</td><td>0.847</td><td>0.867</td><td>0.863</td><td>0.911</td></tr><tr><td>SD2_delta, CorDim_delta, SD2_S2</td><td>Neuronal network</td><td>0.147</td><td>0.454</td><td>1</td><td>0.5</td><td>0.847</td><td>0.867</td><td>0.887</td><td>0.911</td></tr><tr><td>sd_HR_delta, SD2_S2, mean_HR_S2</td><td>Neuronal network</td><td>0.208</td><td>0.77</td><td>1</td><td>0.521</td><td>0.825</td><td>0.869</td><td>0.875</td><td>0.901</td></tr><tr><td>SD2_delta, sd_HR_delta, mean_HR_S2</td><td>GBM</td><td>0.112</td><td>0.331</td><td>0.95</td><td>0.563</td><td>0.881</td><td>0.869</td><td>0.878</td><td>0.9</td></tr><tr><td>SD2_delta, sd_HR_delta, CorDim_delta</td><td>Neuronal network</td><td>0.162</td><td>0.541</td><td>1</td><td>0.375</td><td>0.806</td><td>0.822</td><td>0.833</td><td>0.884</td></tr><tr><td>sd_HR_delta, CorDim_delta, SD2_S2</td><td>Neuronal network</td><td>0.212</td><td>0.748</td><td>1</td><td>0.417</td><td>0.781</td><td>0.827</td><td>0.834</td><td>0.871</td></tr><tr><td colspan="10"><strong>4 variables</strong></td></tr><tr><td>SD2_delta, CorDim_delta, SD2_S2, mean_HR_S2</td><td>GBM</td><td>0.076</td><td>0.245</td><td>0.96</td><td>0.813</td><td>0.935</td><td>0.933</td><td>0.917</td><td>0.947</td></tr><tr><td>sd_HR_delta, CorDim_delta, SD2_S2, mean_HR_S2</td><td>Neuronal network</td><td>0.23</td><td>0.83</td><td>1</td><td>0.625</td><td>0.898</td><td>0.911</td><td>0.903</td><td>0.942</td></tr><tr><td>SD2_delta, sd_HR_delta, SD2_S2, mean_HR_S2</td><td>GBM</td><td>0.097</td><td>0.302</td><td>0.95</td><td>0.813</td><td>0.942</td><td>0.936</td><td>0.935</td><td>0.94</td></tr><tr><td>SD2_delta, sd_HR_delta, CorDim_delta, mean_HR_S2</td><td>GBM</td><td>0.115</td><td>0.348</td><td>0.971</td><td>0.771</td><td>0.885</td><td>0.891</td><td>0.913</td><td>0.916</td></tr><tr><td>SD2_delta, sd_HR_delta, CorDim_delta, SD2_S2</td><td>Neuronal network</td><td>0.213</td><td>0.708</td><td>1</td><td>0.5</td><td>0.847</td><td>0.867</td><td>0.863</td><td>0.911</td></tr><tr><td colspan="10"><strong>5 variables</strong></td></tr><tr><td>SD2_delta, sd_HR_delta, CorDim_delta, SD2_S2, mean_HR_S2</td><td>GBM</td><td>0.092</td><td>0.279</td><td>1</td><td>0.708</td><td>0.878</td><td>0.913</td><td>0.909</td><td>0.931</td></tr></tbody></table>

we found that orthostatism was paradoxically associated with significantly lower low frequency spectral power (LF_power as explicated by the main effects interaction) while high frequency spectral power (HF_power as explicated by the main effects interaction) tends to be lower. Orthostatism is known to physiologically induce a temporary adaptative response increasing sympathetic and decreasing parasympathetic activities. Oppositely to previous descriptions⁹,²⁹, our data do not reveal any major alteration in the sympathovagal response, as expressed by LF_prc, LF_nu, HF_prc, HF_nu and LF_HF ratio, to orthostatism. These normal responses are explicitly evident in Table 3 of the results. LF_power_nu and even more so HF_power, HF_power_prc, HF_power_nu, as well as LF_HF_power, are significantly modified during orthostatism compared to decubitus (position effect, p < 0.001). No clinical form effect (group effect) was observed but only main effects interactions (main effects interactions, p < 0.001). The effect of orthostatism was different according to clinical form. Yet, indexes from RRI time series spectral analysis are known to be markers of the autonomic nervous system influences on the
sinus node rhythmic activity. Particularly LF and HF are thought to reflect sympathovagal balance as well as sympathetic and parasympathetic sinus influences, respectively, when expressed in normalized units. According to first HRV guidelines¹³, one could interpret our results as a sympathovagal balance functionality. On the contrary, there are at least two factors that should prompt us to approach the interpretation more judiciously. Firstly, a recent and highly regarded review³⁰ emphasizes that scientific evidence in the literature supports a predominantly parasympathetic-determined Heart Rate Variability (HRV) when assessed through spectral analysis indexes. This implies that, based on the results of HFnu and LFnu, as well as others spectral analysis outcomes, it appears that parasympathetic-determined HRV does not play a major role in partitioning the clinical forms of SSc patients. Secondly, the mathematical constructs of normalized HRV spectral metrics can lead to a confusing interpretation, depending on whether LF and HF relatively vary with LFnu = LF/(TP - VLF) ≈ LF/(LF + HF) and HFnu = HF/(TP - VLF) ≈ HF/(LF + HF). According to these knowledges and our observations, it is difficult to interpret our results as a simple and unique parasympathetic blunted response (withdrawal) to transition from decubitus to orthostatic conditions. It seems that in orthostatism, sympathetic (LF_power (ms2), LF_power_prc (%)) modulation of sinus activity is lower in dcSSc compared to lcSSc that is in favor of an unpaired adaptation to orthostatism all the more so as it should increase. As LFnu = LF/(TP - VLF), this impairment seems to be supported by modifications impacting LF, TP, but also VLF that is said to reflect mid- and long-term renin-angiotensin system, blood pressure regulation, thermoregulation, and others long-term humoral and hormonal activities¹³. And second, Poincaré plot short term variability index SD1 was not significantly different between dcSSc and lcSSc groups (no clinical form effect) but orthostatism had differential effects according the groups (main effects interactions, p < 0.002). SD1, that is the standard deviation of the plot data on the minor axis is mathematically equal to root mean square of successive differences of RRI duration described the instantaneous beat-to-beat variability that is neurally supported by the only vagal nerve³¹. Oppositely, Poincaré plot long term variability index SD2 that is said to reflect sympathetic and baroreflex HR modulation³²,³³ was significantly different between dcSSc and lcSSc groups during orthostatism but decubitus (main effects interactions, p < 0.034). All these elements converge with the exclusion of a simple vagal-mediated impairment differential main effect between lcSSc and dcSSc, but mostly toward a baroreflex loop complex impairment. Probably for this reason, metrics of HRV linear dynamics lacked to clearly distinguish dcSSc and lcSSc subtypes.

## HR nonlinear dynamics and system failure

Three metrics quantifying the nonlinear properties of HR dynamics have been found to be more deeply impacted in dcSSc vs. lcSSc groups by orthostatic stress according to stand-up-induced changes results: SD2, ApEn, and CorDim. SD2 that is the long/major axe of the elliptic fit of Poincaré plot is also said to reflect nonlinear long-term RR modulation³¹,³². ApEn reflects the complex structure of a time series and the complex behavior of the system generating the time series, i.e., the system complexity and its components interactions³⁴⁻³⁸. CorDim quantifies the time self-similarity of a signal, i.e., of the RRI time series. Thus, CorDim can be considered as an estimation of the lower threshold of the number of degrees of freedom for the underlying system generating the observed data and is typically used as an index of the overall complexity of the system dynamics estimated from a time series²¹. In healthy people, stand-up induces a loss of HRV because of the arterial baroreflex, an adaptive autonomic nervous system response associating sympathetic activation and parasympathetic withdrawal³⁹. But the complexity of a physiological system is increased when its components are interacting to adequate physiological response³⁴. Taken together, the combined stand-up-induced changes of these three metrics that explore complementary nonlinear properties highlight that the cardiovascular homeostatic system is impaired in dcSSc but not in lcSSc. Stand-up stimulus should increase nonlinearity of HR dynamics, which is a wellness index of physiological systems. But in our study, stand-up decreased HR nonlinearity and complexity, which is an impairment and failure marker of dynamic systems including physiological complex systems³⁴. Accordingly, our results converge to highlight that orthostatism stress lowers dramatically physiological flexibility (adaptability) in dcSSc patients compared to lcSSc. DcSSc are less adaptive to orthostatism than lcSSc. These results seem to be in accordance with previous data on baroreflex impairment and orthostatism intolerance in SSc patients. Normal vasoconstrictive response to an increase in venous transmural pressure is almost abolished in tissues with sclerodermic changes, probably due to sympathetic neuropathy⁴⁰. Besides, advanced and fibrotic forms of SSc are associated with a blunted orthostatic stress response⁹. Finally, orthostatic stress intolerance is significantly associated with microvascular damage as assessed by videocapillaroscopy⁴¹. Moreover, in addition to the above-mentioned pathophysiological elements, we cannot exclude but have even to retain the potential role of the pulmonary abnormalities mainly observed in the dcSSc patients and characterizing restrictive lung disease in the loss of HRV nonlinear properties. It is well-known that lung expansion magnitude determines respiratory sinus arrhythmia that is a major component of heart rate variability through vagal activity,⁴²⁻⁴⁴. Now, we have precisely documented that dcSSc exhibited more frequently restrictive syndrome than lcSSc and that the presence of a restrictive syndrome was correlated with HRV differences at rest, during orthostatism, and changes induced by stand-up test.

## SSc subtypes and SSc phenotyping

SSc phenotyping into dcSSc and lcSSc remains a daily challenge for practitioners as clinical form phenotype is a determinant of prognosis and care strategy of the disease¹. In 2019, Sobanski et al. used 24 clinical and serological variables to cluster 6927 SSc patients in an elegant data analysis study⁴⁵. The first 2 clusters overlapped the usual dcSSc and lcSSc subtypes but authors showed that these 2 clusters may not capture the complete heterogeneity of the disease, several other clusters being associated with differential survival rates. They concluded that organ damage should be taken into consideration. Among other criteria that could help SSc phenotyping, HRV has been previously studied³⁻⁶,⁹. Nevertheless, no major difference had been previously robustly documented and
no operational tool has been developed and tested for a real-life application. Through this study we provided SSc subtypes clinical forms classifying models based on HRV parameters usable at the patient's bedside. Several models are provided, more or less complex, but reasonably easy to compute for a daily use, with different performances allowing the clinician to choose the most appropriate for his medical purpose. The easiest models are that provided by our secondary outcomes. HRV can model SSc skin extent trough linear regressions and model respiratory function, mainly the presence of a restrictive syndrome, through logistic regressions. The first being part of the very definition of SSc clinical forms, the second being a highly known pathological status associated with dcSSc. Nevertheless, the efficient models ($R^2 > 0.5$) are quickly complex in terms of number of included variables and interpretability. As shown in Table 2, most of the lung function test variables were significatively different in lcSSc and dcSSc. But most of these variables are correlated and mostly were highly predictable by linear regression models (results not shown). Because of this data high collinearity, probably because of a causal link that could be a limited stretch reflex from the restrictive lungs (neither demonstrated nor documented in this work), we had to provide a specific methodology to extract the most informative variables that we performed through features engineering step. To our knowledge, we propose for the first time a whole functional test of cardiovascular function, usable at the patient's bedside that includes: (1) a 15-min (10 min decubitus, stan-up, 5 min orthostatism) 1-lead numerical ECG recording, (2) the estimation of 1 to 5 HRV metrics easy to compute, and (3) a simple classifying model using machine learning approaches. Locally, we have chosen the GBM 4-variables classifying model including SD2_delta, CorDim_delta, SD2_S2, and mean_HR_S2 HRV indexes as it has the best specificity and F1-score compromise. Nevertheless, the 1-variable model including the only CorDim has the best sensibility-specificity combination while its F1-score is the lowest (0.82) of the tested models that is, in data science an honorable performance. Finally, the choice of the model has to be driven by the underlying medical question and expectation.

## Limits

### Physiological interpretation

A first potential limitation of this study is that we did not directly assess sympathetic and parasympathetic tones by neurography. Physiological and physio-pathological interpretation of the models is then not obvious. We used indirect indexes (spectral powers) from spectral analysis of RRI time series that are known to capture only linear properties of HR dynamics. Nevertheless, in stationary conditions, LF_power and HF_power, expressed in ms² as well as in normalized units, are largely thought to reflect sympathetic and parasympathetic modulatory effects on sinus node activity¹³. As previously discussed, compiling knowledge on effects of sympathetic and vagal activity on HRV suggests that the HRV power spectrum, including its low frequency component, is mainly determined by the parasympathetic system³⁰. Thus, LF-power and HF_power are said to be of limited interest in some physiological conditions, especially during interacting processes that, by definition, lead to nonlinear behavior⁴⁶,⁴⁷. In relation to the complexity of the sinus node activity modulation system, predominantly nonlinear behavior must be assumed, explored, and interpreted.

### Length of time series

A second potential limitation of this study is the number of points constituting the analyzed time series. On one hand, the length of RRI time series used to perform linear time- and frequency-domains HRV analysis is usually set to 5 min¹³. The value of the scaling exponent is mostly defined by the dynamics of the short-time variability. On the other hand, and oppositely, exploring nonlinear dynamics as the dimensionality of the space spanned by the data, particularly when using CorDim, requires long time series to be reliably computed⁴⁸. Time series are supposed to include $10^{CorDim}$ data points, i.e., around 10,000 data points should be used. Using only and approximately 400 data points per time series, we can't consider that CorDim we computed represent robustly and reliably the whole concept and the underlying properties defined by the correlation dimension. But decidedly, the computation we performed (that of CorDim, i.e., the correlation sum) on 400 data points time series led to a metrics that was statistically characterized as sensitive, reliable, stable and discriminant marker. Accordingly, we then highlight that the interpretation of the time behavior of the RRI time series and HR dynamics should be made with caution: The CorDim changes we measured are not necessarily the results of the changes in the scaling behavior of the heart rate dynamics. Specific studies should test this hypothesis.

### Design of the study

A third limitation of the study could be the necessary precaution to take before generalizing these results. Most of the previous studies were designed as sex and age matched groups analysis. On the contrary, we have taken 60 consecutive patients and did not perform specific matching. As shown in Table 1, no age or sex difference appeared but height was significantly different in dcSSc and lcSSc. Three points needs to be addressed to consolidate conclusions: even if our sample size was sufficient for data analysis and to provide, in fine, a validated set of models, sample size has not been calculated on all variables of interest of the study that has not been designed accordingly; no matching has been made on any variable; and no early-stage disease group has been identified since all patient had already their diagnosis. Ormea et al.²⁵ showed that autonomic degeneration was independent of tissue and vascular alterations, as severe nerve damage was already present in apparently healthy skin. He suggested that ANS modification takes place before vascular damage and tissue fibrosis occur. Moreover, Ormea stated that ANS dysfunction is a main factor in the development of the disease. Accordingly, further studies on a larger dataset are needed to confirm and generalize our first results, particularly with the aim to use the classifying models as a predicting tool when used at the early-stage of the disease.
Conclusion, contribution and prospects

In conclusion, this study provides high performance classifying models able to distinguish dcSSc from lcSSc based on only 1 to 5 HRV indexes used as markers of autonomic integrated influences on cardiac activity. The best classifying model are a gradient boosting machine and a neural network reaching a 0.94 F₁-score that include 4 HRV indexes capturing nonlinear properties of HR dynamics. Even if the division of SSc in dcSSc and lcSSc gives mainly a clue of clinical and prognostic stratification, knowing as soon as possible the clinical subtype may impact patient's care. Easy to compute at the patient's bedside, these classifying models need external validation step and have to be tested to predict subtype and clinical form at the early-stage of the disease before an intensive daily medical use.

Data availability

The data that support the findings of this study are available from Assistance Publique-Hôpitaux de Marseille, Marseille, France, but restrictions apply to the availability of these data and so are not publicly available. Data are however available from the authors upon reasonable request and with permission of Assistance Publique-Hôpitaux de Marseille, Marseille, France.

Received: 21 September 2023; Accepted: 24 April 2024

Published online: 15 May 2024

References

1. Volkmann, E. R., Andréasson, K. & Smith, V. Systemic sclerosis. *Lancet* **401**(10373), 304–318 (2023).

2. Hermosillo, A. G., Ortiz, R., Dábague, J., Casanova, J. M. & Martínez-Lavin, M. Autonomic dysfunction in diffuse scleroderma vs CREST: An assessment by computerized heart rate variability. *J. Rheumatol.* **21**(10), 1849–1854 (1994).

3. Morelli, S. et al. Twenty-four-hour heart period variability in systemic sclerosis. *J. Rheumatol.* **23**(4), 643–645 (1996).

4. Ferri, C. et al. Autonomic dysfunction in systemic sclerosis: Time and frequency domain 24 hour heart rate variability analysis. *Br. J. Rheumatol.* **36**(6), 669–676 (1997).

5. Pancera, P. et al. Autonomic nervous system dysfunction in sclerodermic and primary Raynaud's phenomenon. *Clin. Sci.* **96**(1), 49–57 (1999).

6. Colaci, M. et al. Reduction of carotid baroreceptor sensitivity in systemic sclerosis. *Clin. Exp. Rheumatol.* **40**(10), 1964–1969 (2022).

7. Grossman, P. Respiratory sinus arrhythmia (RSA), vagal tone and biobehavioral integration: Beyond parasympathetic function. *Biol. Psychol.* **25**(186), 108739 (2023).

8. Wozniak, J. et al. Evaluation of heart rhythm variability and arrhythmia in children with systemic and localized scleroderma. *J. Rheumatol.* **36**(1), 191–196 (2009).

9. Rodrigues, G. D. et al. Cardiac autonomic modulation at rest and during orthostatic stress among different systemic sclerosis subsets. *Eur. J. Intern. Med.* **66**, 75–80 (2019).

10. Delliaux, S., Delaforge, A., Deharo, J. C. & Chaumet, G. Mental workload alters heart rate variability, lowering non-linear dynamics. *Front. Physiol.* **14**(10), 565 (2019).

11. van den Hoogen, F. et al. 2013 classification criteria for systemic sclerosis: An American college of rheumatology/European league against rheumatism collaborative initiative. *Ann. Rheum. Dis.* **72**(11), 1747–1755 (2013).

12. LeRoy, E. C. et al. Scleroderma (systemic sclerosis): Classification, subsets and pathogenesis. *J. Rheumatol.* **15**(2), 202–205 (1988).

13. Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology. Heart rate variability: Standards of measurement, physiological interpretation and clinical use. *Circulation.* **93**(5), 1043–1065 (1996).

14. Billman, G. E. Heart rate variability—a historical perspective. *Front. Physiol.* **29**(2), 86 (2011).

15. Akselrod, S. et al. Power spectrum analysis of heart rate fluctuation: A quantitative probe of beat-to-beat cardiovascular control. *Science.* **213**(4504), 220–222 (1981).

16. Hoshi, R. A., Pastre, C. M., Vanderlei, L. C. & Godoy, M. F. Poincaré plot indexes of heart rate variability: Relationships with other nonlinear variables. *Auton. Neurosci.* **177**(2), 271–274 (2013).

17. Mourot, L. et al. Decrease in heart rate variability with overtraining: Assessment by the Poincaré plot analysis. *Clin. Physiol. Funct. Imaging.* **24**(1), 10–18 (2004).

18. De Vito, G., Galloway, S. D., Nimmo, M. A., Maas, P. & McMurray, J. J. Effects of central sympathetic inhibition on heart rate variability during steady-state exercise in healthy humans. *Clin. Physiol. Funct. Imaging.* **22**(1), 32–38 (2002).

19. Acharya, R. U., Lim, C. M. & Joseph, P. Heart rate variability analysis using correlation dimension and detrended fluctuation analysis. *ITBM-RBM* **23**, 333–339 (2002).

20. Skinner, J. E., Zebrowski, J. J. & Kowalik, Z. J. New nonlinear algorithms for analysis of heart rate variability: Low dimensional chaos predicts lethal arrhythmias. In *Nonlinear Analysis of Physiological Data* (eds Kantz, H. et al.) 129–166 (Springer, 1998).

21. Sammer, G. Heart period variability and respiratory changes associated with physical and mental load: Non-linear analysis. *Ergonomics.* **41**(5), 746–755 (1998).

22. Grassberger, P. & Procaccia, I. Characterization of strange attractors. *Phys. Rev. Lett.* **50**, 346–349 (1983).

23. Pons, J. F. et al. Heart rhythm characterization through induced physiological variables. *Sci. Rep.* **7**(1), 5059 (2017).

24. Michel, P., Ngo, N., Pons, J. F., Delliaux, S. & Giorgi, R. A filter approach for feature selection in classification: Application to automatic atrial fibrillation detection in electrocardiogram recordings. *BMC Med. Inform. Decis. Mak.* **21** (Suppl 4), 130 (2021).

25. Ormea, F. Studio comparativo del sistema neurovegetativo periferico nella sclerodermia diffusa e nel tessuto connettivale alterato di alcune dermatosi. *Dermatologia* **105**, 8–17 (1952).

26. deBoer, R. W., Karemaker, J. M. & Strackee, J. Hemodynamic fluctuations and baroreflex sensitivity in humans: A beat-to-beat model. *Am. J. Physiol.* **253**(3 Pt 2), H680–H689 (1987).

27. Zlatanovic, M. et al. Cardiac mechanics and heart rate variability in patients with systemic sclerosis: The association that we should not miss. *Rheumatol. Int.* **37**(1), 49–57 (2017).

28. Tadic, M. et al. Systemic sclerosis impacts right heart and cardiac autonomic nervous system. *J. Clin. Ultrasound.* **46**(3), 188–194 (2018).

29. Rodrigues, G. D. et al. Sympatho-vagal dysfunction in systemic sclerosis: A follow-up study. *Life* **13**(1), 34 (2022).

30. Reyes del Paso, G. A., Langewitz, W., Mulder, L. J., van Roon, A. & Duschek, S. The utility of low frequency heart rate variability as an index of sympathetic cardiac tone: A review with emphasis on a reanalysis of previous studies. *Psychophysiology.* **50**(5), 477–487 (2013).

31. Tulppo, M. P., Mäkikallio, T. H., Takala, T. E., Seppänen, T. & Huikuri, H. V. Quantitative beat-to-beat analysis of heart rate dynamics during exercise. *Am. J. Physiol.* **271**(1 Pt 2), H244–H252 (1996).

32. Huikuri, H. V. et al. Abnormalities in beat-to-beat dynamics of heart rate before the spontaneous onset of life-threatening ventricular tachyarrhythmias in patients with prior myocardial infarction. *Circulation.* **93**(10), 1836–1844 (1996).
33. Guzik, P. et al. Correlations between the Poincaré plot and conventional heart rate variability parameters assessed during paced breathing. *J. Physiol. Sci.* **57**(1), 63–71 (2007).

34. Baranger, M. *Complexity, Chaos, and Entropy* (New England Complex Systems Institute, 2000).

35. Pincus, S. M. Approximate entropy as a measure of system complexity. *Proc. Natl. Acad. Sci. USA.* **88**(6), 2297–2301 (1991).

36. Pincus, S. M., Gladstone, I. M. & Ehrenkranz, R. A. A regularity statistic for medical data analysis. *J. Clin. Monit.* **7**(4), 335–345 (1991).

37. Pincus, S. M. & Goldberger, A. L. Physiological time-series analysis: What does regularity quantify?. *Am. J. Physiol.* **266**(4 Pt 2), H1643–H1656 (1994).

38. Pincus, S. M. Greater signal regularity may indicate increased system isolation. *Math Biosci.* **122**(2), 161–181 (1994).

39. Butler, G. C., Yamamoto, Y., Xing, H. C., Northey, D. R. & Hughson, R. L. Heart rate variability and fractal dimension during orthostatic challenges. *J. Appl. Physiol.* **75**(6), 2602–2612 (1993).

40. Henriksen, O., Kristensen, J. K. & Wadskov, S. Local regulation of blood flow in subcutaneous tissue in generalized scleroderma. *J. Investig. Dermatol.* **68**(5), 318–321 (1977).

41. Masini, F. et al. Autonomic nervous system dysfunction correlates with microvascular damage in systemic sclerosis patients. *J. Scleroderma Relat. Disord.* **6**(3), 256–263 (2021).

42. Saul, J. P., Berger, R. D., Chen, M. H. & Cohen, R. J. Transfer function analysis of autonomic regulation. II. Respiratory sinus arrhythmia. *Am. J. Physiol. Heart Circ. Physiol.* **256**(1), H153–H161 (1989).

43. Grossman, P., Wilhelm, F. H. & Spoerle, M. Respiratory sinus arrhythmia, cardiac vagal control, and daily activity. *Am. J. Physiol.-Heart Circ. Physiol.* **287**(2), H728–H734 (2004).

44. Ritz, T. Studying noninvasive indices of vagal control: The need for respiratory control and the problem of target specificity. *Biol. Psychol.* **80**, 158–168 (2009).

45. Sobanski, V. et al. Phenotypes determined by cluster analysis and their survival in the prospective European scleroderma trials and research cohort of patients with systemic sclerosis. *Arthritis Rheumatol.* **71**(9), 1553–1570 (2019).

46. Castiglioni, P. & Parati, G. Present trends and future directions in the analysis of cardiovascular variability. *J. Hypertens.* **29**(7), 1285–1288 (2011).

47. Silva, L. E. V. et al. Comparison between spectral analysis and symbolic dynamics for heart rate variability analysis in the rat. *Sci. Rep.* **7**(1), 8428 (2017).

48. Säkki, M., Kalda, J., Vainu, M. & Laan, M. What does measure the scaling exponent of the correlation sum in the case of human heart rate?. *Chaos* **14**(1), 138–144 (2004).

## Author contributions

SD conceptualized and designed the study, performed the HRV analysis and exploratory data processing, analyzed the results, wrote the paper and led the research; AS aggregated medical records, managed the clinical data acquirement and build the medical database; AE performed all data engineering from descriptive statistics to models' performances comparison, analyzed and discussed the results, and partially wrote the methodological part of the paper; AB discussed the results and commented on the manuscript; QGP discussed the results and commented on the manuscript; FB contributed to patient's recruitment and to the organization of the research, discussed the results and commented on the manuscript; BG formulated the medical need, managed the clinical data acquirement and build the medical database, discussed the results and commented on the manuscript.

## Funding

The project leading to this publication has received funding via Läennec institute from the Excellence Initiative of Aix-Marseille Université-A*Midex, a French "Investissements d'Avenir programme" AMX-21-IET-017.

## Competing interests

The authors declare no competing interests.

## Additional information

Correspondence and requests for materials should be addressed to S.D.

Reprints and permissions information is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_6.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024