Clusters of Survivors of COVID-19 Associated Acute Respiratory Failure According to Response to Exercise

COVID-19 survivors are associated with acute respiratory failure (ARF) and show a high prevalence of impairment in physical performance. The present studied aimed to assess whether we may cluster these individuals according to an exercise test. The presented study is a retrospective analysis of 154 survivors who were admitted to two hospitals of Istituti Clinici Scientifici Maugeri network, Italy. Clinical characteristics, walked distance, heart rate (HR), pulse oximetry (SpO2), dyspnoea, and leg fatigue (Borg scale: Borg-D and Borg-F, respectively) while performing the six-minute walking test (6MWT) were entered into unsupervised clustering analysis. Multivariate linear regression identified variables that were informative for the set of variables used for cluster definition. Cluster 1 (C1: 86.4% of participants) and Cluster 2 (C2: 13.6%) were identified. Compared to C1, the individuals in C2 were significantly older, showed significantly higher increase in fatigue and in dyspnoea, greater reduction in SpO2, and a lower HRpeak during the test. The need of walking aids, time from admission to acute care hospitals, age, body mass index, endotracheal intubation, baseline HR and baseline Borg-D, and exercise-induced SpO2 change were significantly associated with the variables that were used for cluster definition. Different characteristics and physiological parameters during the 6MWT characterise survivors of COVID-19-associated ARF. These results may help in the management of the long-term effects of the disease.


Introduction
The SARS-coronavirus-2 disease 19 (SARS-CoV-2) pandemic has had dramatic effects throughout the world, with more than two hundred, thirty million people infected and more than four and half million casualties worldwide, requiring comprehensive management strategies in the acute as well as in the post-acute phases of the pandemic [1,2]. In addition to the consequences on lung function [3], a high prevalence of impairment in physical performance is reported in survivors, who may suffer from fatigue and/or muscle weakness, exercise-induced dyspnoea, sleep difficulties, and anxiety and/or depression up to 6 months after infection [4][5][6][7]. In addition, persistent exercise-induced desaturation (EID) was reported in up to 43% of survivors of acute respiratory failure (ARF) due to COVID-19 pneumonia [8], and these individuals, when normoxaemic at rest with EID, may show Int. J. Environ. Res. Public Health 2021, 18, 11868 2 of 11 alterations in lung function, exercise capacity, and symptoms similar to individuals with interstitial lung diseases but that are more severe than those of individuals with chronic obstructive pulmonary disease (COPD) and EID [9].
We wondered whether these individuals might be characterized on the basis of physiological response to symptoms during an exercise test. Cluster analysis consists of methods to identify subgroups of data. By assigning items into the subgroups that are identified, clustering can help us to understand the characteristics of the patterns present in data. These techniques are commonly used in a wide range of areas, including medicine, psychology, market research, and bioinformatics [10].
Therefore, the aim of this study was to assess whether survivors of COVID-19 with associated ARF, and who were normoxaemic at rest, could be clustered according to exercise capacity and physiological and symptom response to an exercise test.

Materials and Methods
This retrospective case series study was approved by the Istituti Clinici Scientifici (ICS) Maugeri, IRCCS Ethics Committee (2440 CE). As a retrospective study, the individuals did not need to provide any specific written informed consent; however, at admission they gave informed consent for the scientific use of their clinical data in advance. As this was a retrospective analysis, the study was not registered.

Participants
The participants included individuals who were consecutively admitted between 15th April and 30th May 2020 to the Respiratory Units of ICS Maugeri, IRCCS hospitals of Lumezzane (Brescia) and Pavia-Italy, referral institutions for pulmonary rehabilitation, diagnosis, and care of post-acute and chronic conditions [11,12]. These individuals had been transferred from acute care hospitals where they had undergone invasive or non-invasive (NIV) ventilation, high flow oxygen therapy, and/or oxygen supplementation [13].
The inclusion criteria were negative Real Time PCR (RT-PCR) test for SARS-CoV-2 at the time of evaluation, ability to perform the six-minute walking test (6MWT) at discharge from hospital, and to be normoxaemic at rest.
The exclusion criteria were the need of oxygen supplementation at the time of evaluation or the need of Long-Term Oxygen Therapy (LTOT), the presence of neurological or orthopaedic conditions (chronic or new onset), and inability or unwillingness to undergo evaluations.

Measurements
The following data were recorded from the ICS clinical data repository: demographics, anthropometrics, comorbidities, length of stay (LoS) in acute care and in our hospitals, history of endotracheal intubation (EI) and mechanical ventilation, NIV or oxygen supply in the acute care hospitals, presence of tracheostomy, need of walking aids, Motor Barthel index [14], and serological parameters (C-reactive protein [CRP], D-dimer, ferritin, albumin, haemoglobin) at admission.
With the appropriate safety procedures and while wearing appropriate personal protective equipment [15], the participants underwent arterial blood gases under room air at admission.
Exercise tolerance was assessed by means of the 6MWT according to accepted standards [16]. The 6MWT is the most worldwide used and accepted field test to assess exercise capacity in different conditions and setting including post-COVID-19 and during rehab [5]. The walked distance was expressed as meters and as the percent of the reference values [17]. At the beginning and at the end of the walking cycle, subjective sensations of dyspnoea (Borg-D) and leg fatigue (Borg-F) were assessed by means of the modified Borg scale [18]. Peripheral oxygenation (SpO 2 ) and heart rate (HR) were monitored using a pulse oximeter (VintusWalk, VYAIRE MEDICAL, INC, NY, USA). Baseline and SpO 2 nadir and baseline and HR peak were also recorded. Exercise-induced desaturation was defined as SpO 2 nadir -baseline SpO 2 > 4% during the 6MWT [19]. Heart rate and % of maximal predicted (220-age) was also recorded.

Statistical Analysis
Continuous and discrete numeric variable distribution was described as the median and interquartile range (25th, 75th percentiles, IQR) since most of them deviated from the normal distribution based on visual inspection of histograms. Categorical variable distribution was described as the absolute and relative frequency (%). Change in continuous and discrete numeric variables between baseline and end of the 6MWT were calculated as end-exercise-baseline values.
Unsupervised hierarchical agglomerative clustering analysis was performed using the Euclidean distance between the participants and the complete linkage function after data centring and scaling using the mean and standard deviation (SD) of each variable's distribution. The information for the distance during the 6MWT, % predicted, HR change (HR peak -baseline HR), SpO 2 change (SpO 2 nadir -baseline SpO 2 ), Borg-D change (end-exercise Borg-D-baseline Borg-D), and Borg-F change (end-exercise Borg-F-baseline Borg-F) were used for cluster identification on the overall dataset and included participants from both study centres. The optimal number of clusters was identified as the one guaranteeing the highest silhouette coefficient value by testing different numbers of clusters (from 1 to 20).
The Spearman correlation coefficient r and corresponding 95% confidence interval (95% CI) were calculated to estimate the strength of the correlation between continuous and discrete numeric variables. The Wilcoxon rank sum test was applied to compare numeric variable distribution between groups. Pearson chi-square test with p-values computation using the Monte Carlo simulation (10,000 simulations) was applied to test for association between categorical variables.
Multivariate linear regression with stepwise forward features selection based on the Bayesian Information Criterion (BIC) was applied to identify the most informative subset of variables associated with the response variable of interest. Before multivariate linear regression fitting, missing values of numeric continuous and discrete variables with a missing data fraction ≤20% were imputed by the median value of the corresponding distribution and the missing values of the categorical variables by the most frequent value of the corresponding frequency distribution. Variables that were characterized by a missing data fraction >20% were not included in the analysis. A starting set of variables including the hospital centre and the baseline value of each response variable was forced to be included into the model. The set of potentially informative variables based on clinical knowledge was represented by LoS in acute hospitals (days), LoS in study centres (days), time from acute hospital admission (days), age (years), body mass index (BMI, Kg/m 2 ), arterial oxygen tension (PaO 2 , mmHg), EI in acute hospitals (yes, no), NIV in acute hospitals (yes, no), SpO 2 change (%), walking aids requirement (yes, no), and CRP (mg/dL). All analyses were performed using the R statistical software tool (www.r-project.org; version 4.0.5), clustering analyses were performed by the NbClust function implemented in the R package called NbClust, and 95% CIs of the Spearman correlation coefficient r were computed by the spearman.ci function implemented in the RVAideMemoire package. Functions corresponding to the remaining statistical tests used in the analyses were implemented in the stats package.

Patients' Characteristics
Out of the 196 individuals who were admitted during the study period, the data from a total number of 154 participants with complete variable data that could be used for clustering (see Statistical Analysis) were included in the analysis. Demographic and clinical data are shown in Table 1. According to the inclusion criteria, all of participants were normoxaemic at rest. Compared to Pavia, the participants from Lumezzane were more likely to be female, obese, and suffer from congestive heart failure, with a longer LoS in acute hospitals and in study centres and a longer time from admission to acute care hospitals, with a higher Barthel index, CRP, albumin, PaCO 2 , baseline Borg-D, and baseline Borg-F as well as lower D-dimer and ferritin and a lower walked distance at 6MWT and baseline SpO 2 .

Correlation between Clustering Parameters
As expected, all of participants showed an increase in their HR during the 6MWT, with 81 out of 154 (52.2%) showing EID. Out of 154, 129 (83.8%) and 93 (60.4%) participants showed an increase of at least 0.5 points in Borg-D and Borg-F, respectively at the endexercise point.
Measurements of the participants belonging to the different clusters are highlighted with distinct colours and shapes in the scatterplots (grey circles = measures of patients belonging to C1, black crosses = measures of patients belonging to C2). The number reported within each section is the Spearman correlation coefficient describing the strength of the correlation between the variables. As an example, a positive correlation was observed between Borg-D change and Borg-F change (Spearman r = 0.60); patients belonging to C2 tended to be characterized by higher Borg-D change and Borg-F change values compared to the patients belonging to C1.

Clustering Analysis
Through unsupervised analysis it was possible to identify two clusters: Cluster 1 (C1) included 133 participants (86.4%), and Cluster 2 (C2) included the remaining 21 (13.6%). Compared to C1, the individuals in C2 showed a larger increase in Borg-F and Borg-D and a greater reduction in SpO 2 ( Figure 2 and Table 2).

Clustering Analysis
Through unsupervised analysis it was possible to identify two clusters: Cluster 1 (C1) included 133 participants (86.4%), and Cluster 2 (C2) included the remaining 21 (13.6%). Compared to C1, the individuals in C2 showed a larger increase in Borg-F and Borg-D and a greater reduction in SpO2 (Figure 2 and Table 2). From top to bottom, each boxplot represented: the lower non-outlier limit, the 25th percentile, the 50th percentile, the 75th percentile, and the upper non-outlier limit. Each point outside of the non-outlier range represents an outlier with respect to the variables' distribution. p = p-value deriving from the Wilcoxon rank-sum test.
When focusing on the variables that were not used for cluster definition, we observed that compared to C1, the individuals in C2 were older, more likely to be affected by renal failure and COPD, and had higher baseline Borg-D, end-exercise Borg-D, end-exercise Borg-F, and lower HRpeak ( Table 2).  From top to bottom, each boxplot represented: the lower non-outlier limit, the 25th percentile, the 50th percentile, the 75th percentile, and the upper non-outlier limit. Each point outside of the non-outlier range represents an outlier with respect to the variables' distribution. p = p-value deriving from the Wilcoxon rank-sum test.
When focusing on the variables that were not used for cluster definition, we observed that compared to C1, the individuals in C2 were older, more likely to be affected by renal failure and COPD, and had higher baseline Borg-D, end-exercise Borg-D, end-exercise Borg-F, and lower HR peak ( Table 2).

Factors Associated to Variables Used for Clusters Definition
The results of the multivariate feature selection process aimed to identify the most informative subset of factors associated with the set of variables used for cluster definition, which are shown in Table 3 and can be briefly described as follows: (a) The need of walking aids was associated with lower 6MWT and % predicted. Participants from Pavia were more likely to be characterized by higher 6MWT and % predicted. (b) History of EI was associated with a positive shift in terms of HR change distribution, while SpO 2 change, baseline HR, and time from admission to acute care hospitals were associated with a negative shift in terms of HR change distribution. The participants from Pavia were also characterized by a positive shift in terms of HR change distribution. (c) Body mass index was associated with a positive shift in terms of SpO 2 change distribution, and the participants from Pavia centre were also likely to be characterized by a positive shift in terms of SpO 2 change distribution. (d) Baseline Borg-D and the need of walking aids were associated with a positive shift in Borg-D change distribution. SpO 2 change during the test was associated with a negative shift in Borg-D change distribution. (e) The need of walking aids and age were associated with a positive shift in terms of Borg-F change distribution.

Discussion
Survivors of ARF caused by to COVID-19 pneumonia can be distinguished by two clusters that are characterised by different exercise capacities and exercise-induced symp-toms. More than 80% of the participants were included in C1, which was characterised by male participants who were able to achieve a better performance during the 6MWT and higher HR peak . Older participants with more severe symptoms, greater SpO 2 change during 6MWT, and a greater prevalence of renal failure and COPD were included into C2. The need of walking aids, age, BMI, hospital centre, SpO 2 change, baseline HR, time from admission to acute hospital, history of EI, baseline SpO 2 , and baseline Borg-D were informative with respect to 6MWT, % predicted, changes in SpO 2 , changes in HR, and changes by end exercise-induced severity of dyspnoea and fatigue, respectively.
By cluster analysis, our study has shown the ability to distinguish between the different subgroups COVID-19 pneumonia survivors. Cluster analysis consists of exploratory techniques that are widely used in different research fields, including biology and medicine [10,20]. Clustering allows the discovery subgroups of subjects who share similar characteristics. While our study provides interesting information on exercise capacity in survivors of COVID-19-associated pneumonia, our study is also relevant, as it suggests a modality that can be used to analyse the characteristics of these patients in order to add to the panel of common evaluations. In addition, the characterisation of subgroups according to demographic, clinical, physiological characteristics and response to exercise may have potential clinical usefulness. Different characteristics and responses to a simple and easy to perform exercise test would allow the tailoring of comprehensive management strategies, including exercise training, and may also suggest different hypothetical recovery times after an acute hospital admission or different prognostic consequences in this new population, which is now called post-COVID patients. Prospective studies might confirm the evaluation of these perspectives.
Confirming a previous study [9], our participants showed a wide range of exercise capacities, as assessed by 6MWT values ranging from about 50% to 90%, with median values that were less than 70% of those predicted. To what extent the reduction in exercise capacity was due to the infection and/or to clinical consequences such as prolonged hospital LoS and/or mechanical ventilation is still to be elucidated. The participants in the current study were assessed at a median time of 38 days (min = 8, max = 120 days) from acute hospital admission. However, the time from this event was correlated to change in HR during exercise but not to the distance walked (see Table 3).
In the study by Huang et al. [6], in most individuals who recovered from severe COVID-19, dyspnoea scores and exercise capacity improved over time, but in that study [6], individuals who required EI and mechanical ventilation were excluded due to the potential for the consequences of mechanical ventilation itself influencing the factors that were under investigation. In our study, the need for EI was informative of the changes that take place in terms of HR during exercise but were not indicative of exercise capacity itself (see Table 3).
Participants from the two centres were different in terms of some demographic, anthropometric, clinical, and physiological characteristics as well as in terms of responses to the exercise test. However, the contributions of the two centres to the clusters did not differ significantly. Of the variables that were characterized by unbalanced distributions between centres, only the baseline Borg-D was differentially distributed between clusters; thus, it is likely that these variables did not heavily influence cluster definition.

Limitations
As this is a retrospective study, the present work may contain flaws that are associated with such studies, including missing data. Indeed, the pandemic has led to an increase in retrospective publications that have a high level of restrictions [21]. However, we could not wait for well-designed prospective randomized controlled trials before starting interventions in daily clinical practice. Therefore, we feel that with all of the in-built limitations of observational, uncontrolled, or retrospective studies, the scientific community has had to answer the emerging questions posed by the pandemic, which also applied to the field of rehabilitation. Because of this, making timely use of the available data is integral.
We did not study a control population; however, we are confident that the comparison with predicted 6MWD values allowed us to conduct an unbiased evaluation of the results.
The lack of lung function assessment may be criticized. This was an unavoidable consequence of the preventive protection measures imposed by the pandemic. Indeed, the study involved patients who were admitted during the biggest outbreak of the pandemic (April-May 2020), when all activities that could potentially spread the infection were suspended, including lung function studies.
The reported occurrence of commodities in our study must not be considered as a true indication of prevalence, as our patients did not undergo any specific diagnostic tests.

Conclusions
Despite the limitations of a retrospective design, using cluster analysis, our study allowed us to identify subpopulations of survivors of COVID-19-associated pneumonia. These results might help in the management of the long-term effects of the disease but should be confirmed by prospective, controlled studies.

Informed Consent Statement:
Written informed consent has been obtained from the patient(s) to publish this paper.

Data Availability Statement:
Anonymized data materials will be made publicly available at https: //www.zenodo.org/.