Investigating phenotypes of pulmonary COVID-19 recovery: A longitudinal observational prospective multicenter trial

Abstract

Background: The optimal procedures to prevent, identify, monitor, and treat long-term pulmonary sequelae of COVID-19 are elusive. Here, we characterized the kinetics of respiratory and symptom recovery following COVID-19.

Methods: We conducted a longitudinal, multicenter observational study in ambulatory and hospitalized COVID-19 patients recruited in early 2020 (n = 145). Pulmonary computed tomography (CT) and lung function (LF) readouts, symptom prevalence, and clinical and laboratory parameters were collected during acute COVID-19 and at 60, 100, and 180 days follow-up visits. Recovery kinetics and risk factors were investigated by logistic regression. Classification of clinical features and participants was accomplished by unsupervised and semi-supervised multiparameter clustering and machine learning.

Results: At the 6-month follow-up, 49% of participants reported persistent symptoms. The frequency of structural lung CT abnormalities ranged from 18% in the mild outpatient cases to 76% in the intensive care unit (ICU) convalescents. Prevalence of impaired LF ranged from 14% in the mild outpatient cases to 50% in the ICU survivors. Incomplete radiological lung recovery was associated with increased anti-S1/S2 antibody titer, IL-6, and CRP levels at the early follow-up. We demonstrated that the risk of perturbed pulmonary recovery could be robustly estimated at early follow-up by clustering and machine learning classifiers employing solely non-CT and non-LF parameters.

Conclusions: The severity of acute COVID-19 and protracted systemic inflammation is strongly linked to persistent structural and functional lung abnormality. Automated screening of multiparameter health record data may assist in the prediction of incomplete pulmonary recovery and optimize COVID-19 follow-up management.

Funding: The State of Tyrol (GZ 71934), Boehringer Ingelheim/Investigator initiated study (IIS 1199-0424).

Clinical trial number: ClinicalTrials.gov: NCT04416100.

Keywords: COVID-19; computed tomography; epidemiology; global health; human; long COVID; machine learning; medicine; post-COVID-19 syndrome; pulmonary recovery.

Conflict of interest statement

TS, SS, AB, AP, AL, CS, MN, PG, KK, SK, MA, BP, AE, GH, EW, GW, GW, IT, JL No competing interests declared, PT owns his own business, Data Analytics as a Service Tirol, for which he performs freelance data science work. Has also received an honorarium for the study data management, curation and analysis and minor manuscript work. The author has no other competing interests to declare

© 2022, Sonnweber et al.

Figures

Figure 1.. Study inclusion flow diagram and analysis scheme.

Figure 2.. Kinetic of recovery from COVID-19 symptoms.

Recovery from any COVID-19 symptoms was investigated by mixed-effect logistic modeling (random effect: individual; fixed effect: time). Significance was determined by the likelihood ratio test corrected for multiple testing with the Benjamini–Hochberg method, and p-values and the numbers of complete observations are indicated in the plots. (A) Frequencies of individuals with any symptoms in the study cohort stratified by acute COVID-19 severity. (B) Frequencies of participants with particular symptoms. imp.: impaired.

Figure 3.. Kinetic of pulmonary recovery.

Recovery from any lung computed tomography (CT) abnormalities, moderate-to-severe lung CT abnormalities (severity score > 5), and recovery from functional lung impairment were investigated in the participants stratified by acute COVID-19 severity by mixed-effect logistic modeling (random effect: individual; fixed effect: time). Significance was determined by the likelihood ratio test corrected for multiple testing with the Benjamini–Hochberg method. Frequencies of the given abnormality at the indicated time points are presented, and p-values and the numbers of complete observations are indicated in the plots.

Figure 3—figure supplement 1.. Co-occurrence of lung computed tomography (CT) abnormalities, functional lung impairment, and any persistent symptoms.

Numbers and percentages of the study participants with any persistent symptoms, functional lung impairment, or lung CT abnormalities at the consecutive follow-up visits presented in quasi-proportional Venn diagrams. The numbers of participants with CT abnormalities, lung function (LF) impairment, and persistent symptoms are indicated in the diagrams, and the numbers of complete observations are shown under the plots.

Figure 3—figure supplement 2.. Co-occurrence of moderate-to-severe lung computed tomography (CT) abnormalities, functional lung impairment, and any persistent symptoms.

Numbers and percentages of the study participants with any persistent symptoms, functional lung impairment, or moderate-to-severe lung CT abnormalities (severity score > 5) at the consecutive follow-up visits presented in quasi-proportional Venn diagrams. The numbers of participants with CT abnormalities, lung function (LF) impairment, and persistent symptoms are indicated in the diagrams, and the numbers of complete observations are shown under the plots.

Figure 3—figure supplement 3.. Frequency of mild and moderate-to-severe lung computed tomography (CT) abnormalities.

Prognostic value of functional lung impairment and persistent symptoms for prediction of radiological lung abnormalities. (A) Relevance of functional lung impairment and persistent COVID-19 symptoms at predicting any lung CT abnormalities and moderate-to-severe lung CT abnormalities (severity score > 5) at the consecutive follow-up visits. The concordance of the outcome variables was determined by Cohen’s κ coefficient. Statistical significance (κ / = 0) was assessed by two-tailed t-test corrected for multiple testing with the Benjamini–Hochberg method. Kappa with 95% confidence intervals and p values are presented as a heat map. The number of complete observations is indicated in the plot. (B) Percentages of mild (severity score ≤ 5) and moderate-to-severe lung CT abnormalities at the consecutive follow-up visits in the study participants stratified by the severity of acute COVID-19. Statistical significance of frequency differences was determined by χ2 test for trend corrected for multiple testing with the Benjamini–Hochberg method. The number of complete observations is indicated under the plots.

Figure 4.. Risk factors of persistent radiological lung abnormalities.

Association of 52 binary explanatory variables (Appendix 1—table 1) with the presence of any lung computed tomography (CT) abnormalities (A) or moderate-to-severe lung CT abnormalities (severity score > 5) (B) at the 180-day follow-up visit was investigated with a series of univariate logistic models (Appendix 1—table 2). Odds ratio (OR) significance was determined by Wald Z test and corrected for multiple testing with the Benjamini–Hochberg method. ORs with 95% confidence intervals for significant favorable and unfavorable factors are presented in forest plots. Model baseline (ref) and numbers of complete observations are presented in the plot axis text. Q1, Q2, Q3, Q4: first, second, third, and fourth quartile of anti-S1/S2 IgG titer; ICU: intensive care unit.

Figure 5.. Risk factors of persistent functional lung impairment.

Association of 52 binary explanatory variables (Appendix 1—table 1) with the presence of functional lung impairment at the 180-day follow-up visit was investigated with a series of univariate logistic models (Appendix 1—table 2). Odds ratio (OR) significance was determined by Wald Z test and corrected for multiple testing with the Benjamini–Hochberg method. ORs with 95% confidence intervals for the significant favorable and unfavorable factors are presented in a forest plot. Model baseline (ref) and n numbers of complete observations are presented in the plot axis text. Q1, Q2, Q3, Q4: first, second, third, and fourth quartile of anti-S1/S2 IgG titer; CKD: chronic kidney disease.

Figure 6.. Association of incomplete symptom, lung function, and radiological lung recovery with demographic and clinical parameters of acute COVID-19 and early recovery.

Clustering of 52 non-computed tomography (non-CT) and non-lung function binary explanatory variables recorded for acute COVID-19 or at the early 60-day follow-up visit (Appendix 1—table 1) was investigated by partitioning around medoids (PAM) algorithm with simple matching distance (SMD) dissimilarity measure (Figure 6—figure supplement 1, Appendix 1—table 3). The cluster assignment for the outcome variables at the 180-day follow-up visit (persistent symptoms, functional lung impairment, mild lung CT abnormalities [severity score ≤ 5] and moderate-to-severe lung CT abnormalities [severity score > 5]) was predicted by k-nearest neighbor (kNN) label propagation procedure. Numbers of complete observations and numbers of features in the clusters are indicated in (A). (A) Cluster assignment of the outcome variables (diamonds) presented in the plot of principal component (PC) scores. The first two major PCs are displayed. The explanatory variables are visualized as points. Percentages of the data set variance associated with the PC are presented in the plot axes. (B) Five nearest neighbors (lowest SMD) of the outcome variables presented in radial plots. Font size, point radius, and color code for SMD values. Q1, Q2, Q3, Q4: first, second, third, and fourth quartile of anti-S1/S2 IgG titer; GITD: gastrointestinal disease; CKD: chronic kidney disease; ICU: intensive care unit; COPD: chronic obstructive pulmonary disease.

Figure 6—figure supplement 1.. Study feature clustering algorithm.

Clustering of 52 non-computed tomography (non-CT) and non-lung function binary explanatory variables recorded for acute COVID-19 or at the early 60-day follow-up visit (Appendix 1—table 1). (A, B) Comparison of ‘explained’ variances (between-cluster to total sum-of-squares ratio) (A) and cluster stability (mean classification error in 20-fold cross-validation) (B) in clustering of the data set with several algorithms with k = 3 centers/branches (algorithms: K-means; PAM: partitioning around medoids; HCl Ward.D2: hierarchical clustering with Ward.D2 method; distances: SMD: simple matching distance; Jaccard, Dice, and Cosine). (C, D) The optimal number of the feature clusters in clustering with the optimally performing PAM algorithm with SMD dissimilarity measure was determined by the bend of the total within-cluster sum-of-squares curve (C) and confirmed by good stability (low mean classification error) in 20-fold cross-validation (D).

Figure 6—figure supplement 2.. Semi-supervised clustering of mild and moderate-to-severe lung computed tomography (CT) abnormalities, functional lung impairment, and persistent symptoms at the 180-day follow-up with parameters of acute COVID-19 and early convalescence.

Clusters of 52 non-CT and non-lung function binary explanatory variables recorded for acute COVID-19 or at the 60-day follow-up visit (Appendix 1—table 1) were defined by the optimally performing partitioning around medoids (PAM) algorithm and simple matching distance (SMD) dissimilarity measure (Figure 6A, Figure 6—figure supplement 1, Appendix 1—table 3). The cluster assignment for the outcome variables at the 180-day follow-up visit (persistent symptoms, functional lung impairment, mild lung CT abnormalities [severity score ≤ 5], and moderate-to-severe lung CT abnormalities [severity score > 5]) was predicted by k-nearest neighbor (kNN) label propagation procedure. SMD between the features and their cluster assignments are shown in a heat map. The numbers of features in the clusters and the total number of observations are indicated under the plot. CVD: cardiovascular disease; Q1, Q2, Q3, Q4: first, second, third, and fourth quartile of anti-S1/S2 IgG titer; GI: gastrointestinal; PD: pulmonary disease; GITD: gastrointestinal disease; ICU: intensive care unit; COPD: chronic obstructive pulmonary disease; CKD: chronic kidney disease.

Figure 7.. Clustering of the study participants by non-lung function and non-computed tomography (non-CT) clinical features.

Study participants (n = 133 with the complete variable set) were clustered with respect to 52 non-CT and non-lung function binary explanatory variables recorded for acute COVID-19 or at the 60-day follow-up visit (Appendix 1—table 1) using a combined self-organizing map (SOM: simple matching distance) and hierarchical clustering (Ward.D2 method, Euclidean distance) procedure (Figure 7—figure supplement 1). The numbers of participants assigned to low-risk (LR), intermediate-risk (IR), and high-risk (HR) clusters are indicated in (A). (A) Cluster assignment of the study participants in the plot of principal component (PC) scores. The first two major PCs are displayed. Percentages of the data set variance associated with the PC are presented in the plot axes. (B) Presence of the most influential clustering features (Figure 7—figure supplement 2) in the participant clusters presented as a heat map. Cluster #1, #2, and #3 refer to the feature clusters defined in Figure 6. Q1, Q2, Q3, Q4: first, second, third, and fourth quartile of anti-S1/S2 IgG titer; GITD: gastrointestinal disease; CKD: chronic kidney disease; CVD: cardiovascular disease; GI: gastrointestinal; PD: pulmonary disease.

Figure 7—figure supplement 1.. Study participant clustering algorithm.

Clustering of the study participants (n = 133 with the complete variable set) with respect to 52 non-computed tomography (non-CT) and non-lung function binary explanatory variables recorded for acute COVID-19 or at the 60-day follow-up visit (Appendix 1—table 1). The procedure involved clustering of the observations with self-organizing maps (SOM, 4 × 4 hexagonal grid, distances: SMD: simple matching distance, Jaccard, Dice, or Cosine) followed by clustering of the SOM nodes (algorithms: HCl ward.D2: hierarchical clustering with Ward.D2 method; K-means; PAM: partitioning around medoids; distance: Euclidean). Different combinations of observation dissimilarity measures and SOM node clustering algorithms were tested in the search for the optimal clustering algorithm. (A, B) Comparison of ‘explained’ variances (between-cluster to total sum-of-squares ratio) (A) and cluster stability (mean classification error in 20-fold cross-validation) (B) in clustering of the data set with different observation distance measures and SOM node clustering algorithms. (C) Training of the SOM algorithm, mean distance to the winning un as a function of lgorithm iterations is presented. Note the mean distance plateau indicative of the algorithm convergence (D–F) The optimal number of the SOM node clusters in clustering with the optimally performing SOM HCl algorithm with SMD observation dissimilarity measure. The optimal cluster number was determined by the bend of the total within-cluster sum-of-squares curve (D) and confirmed by visual inspection of the HCl dendrogram (E) and good stability (low mean classification error) in 20-fold cross-validation (F).

Figure 7—figure supplement 2.. Impact of specific variables on the quality of participant clustering.

The clusters of participants clusters were defined with the optimally performing self-organizing map (SOM)/HCl algorithm with simple matching distance (SMD) observation dissimilarity measure as presented in Figure 7 and Figure 7—figure supplement 1. The impact of a particular clustering variable was determined by comparing the ‘explained’ clustering variance (between-cluster to total sum-of-squares ratio) between the initial cluster structure and the structure with random resampling of the variable. Differences in the clustering variances for the most influential clustering variables (Δ clustering variance > 0) are presented in the plot. Q1, Q3: first, third quartile of anti-S1/S2 IgG titer; CKD: chronic kidney disease; GI: gastrointestinal; CVD: cardiovascular disease; PD: pulmonary disease; GITD: gastrointestinal disease.

Figure 8.. Frequency of persistent radiological lung abnormalities, functional lung impairment, and symptoms in the participant clusters.

The clusters of study participants were defined by non-lung function and non-computed tomography (non-CT) features as presented in Figure 7. Frequencies of outcome variables at the 180-day follow-up visit (mild [severity score ≤ 5], moderate-to-severe lung CT abnormalities [severity score > 5], functional lung impairment, and persistent symptoms) were compared between the low-risk (LR), intermediate-risk (IR), and high-risk (HR) participant clusters by χ2 test corrected for multiple testing with the Benjamini–Hochberg method. p-Values and numbers of participants assigned to the clusters are indicated in the plots. (A) Frequencies of the outcome features in the participant clusters. (B) Frequencies of specific symptoms in the participant clusters.

Figure 8—figure supplement 1.. Risk of radiological lung abnormalities at the 180-day follow-up in the participant clusters.

The clusters of participants were defined by non-lung function and non-computed tomography (non-CT) clinical features of acute COVID-19 and early convalescence (60-day follow-up visit, Appendix 1—table 1) with the optimally performing HCl algorithm with simple matching distance (SMD) observation dissimilarity measure as presented in Figure 7 and Figure 7—figure supplement 1. (A) Distribution of mild, moderate, severe, and critical acute COVID-19 cases in the participant clusters. Significance of the distribution differences was assessed with χ2 test. The numbers of participants assigned to the clusters are indicated under the plot. (B) Association of the participant cluster assignment (LR: low-risk; IR: intermediate-risk; HR: high-risk cluster) with the risk of any lung CT abnormalities and moderate-to-severe lung CT abnormalities (severity score > 5) at the 180-day follow-up visit was investigated by logistic modeling with and without inclusion of the acute COVID-19 severity effect (severity-adjusted). Odds ratio (OR) significance was determined by Wald Z test and corrected for multiple testing with the Benjamini–Hochberg method. ORs with 95% confidence intervals are presented in forest plots. Numbers of complete observations, outcome events, participants in the clusters, and the acute COVID-19 severity subsets are indicated under the plot.

Figure 9.. Prediction of persistent radiological lung abnormalities, functional lung impairment, and symptoms by machine learning algorithms.

Single machine learning classifiers (C5.0; RF: random forests; SVM-R: support vector machines with radial kernel; NNet: neural network; glmNet: elastic net) and their ensemble (Ens) were trained in the cohort data set with 52 non-computed tomography (non-CT) and non-lung function binary explanatory variables recorded for acute COVID-19 or at the 60-day follow-up visit (Appendix 1—table 1) for predicting outcome variables at the 180-day follow-up visit (any lung CT abnormalities, moderate-to-severe lung CT abnormalities [severity score > 5], functional lung impairment, and persistent symptoms) (Appendix 1—table 4). The prediction accuracy was verified by repeated 20-fold cross-validation (five repeats). Receiver-operating characteristics (ROCs) of the algorithms in the cross-validation are presented: area under the curve (AUC), sensitivity (Sens), and specificity (Spec) (Appendix 1—table 5). The numbers of complete observations and outcome events are indicated under the plots.

Figure 9—figure supplement 1.. Correlation of the machine learning algorithm prediction accuracy.

Machine learning classifiers (C5.0; RF: random forests; SVM-R: support vector machines with radial kernel; NNet: neural network; glmNet: elastic net) were trained in the cohort data set with 52 non-computed tomography (non-CT) and non-lung function binary explanatory variables recorded for acute COVID-19 or at the early 60-day follow-up visit (Appendix 1—table 1) for predicting outcome variables at the 180-day follow-up visit (any lung CT abnormalities, moderate-to-severe lung CT abnormalities [severity score > 5], functional lung impairment, and persistent symptoms) (Figure 9, Appendix 1—table 4). The prediction accuracy was verified by repeated 20-fold cross-validation (five repeats). Pearson’s correlation coefficients of the classifier prediction accuracy in the cross-validation are presented as heat maps. Numbers of complete observations and outcome events are indicated under the plots.

Figure 9—figure supplement 2.. Machine learning model ensembles.

Single machine learning classifiers (C5.0; RF: random forests; SVM-R: support vector machines with radial kernel; NNet: neural network; glmNet: elastic net) were trained as shown in Figure 9. The model ensembles based on the single classifiers were constructed with the glmNet procedure (Appendix 1—table 4). glmNet regression coefficients (β) are presented in the plots. Point and text color correspond to the β value. Numbers of complete observations and outcome events are indicated under the plots.

Figure 9—figure supplement 3.. Prediction of persistent radiological lung abnormalities, functional lung impairment, and symptoms by machine learning algorithms in the training data sets.

Single machine learning classifiers (C5.0; RF: random forests; SVM-R: support vector machines with radial kernel; NNet: neural network; glmNet: elastic net) and their ensembles were trained as shown in Figure 9. Performance of the classifiers in the training data sets was investigated by receiver-operating characteristic (ROC) of the algorithms (AUC: area under the curve; Sens: sensitivity; Spec: specificity, Appendix 1—table 5). Numbers of complete observations and outcome events are indicated under the plots.

Figure 9—figure supplement 4.. Variable importance statistics for prediction of lung computed tomography (CT) abnormalities at the 180-day follow-up by machine learning classifiers.

C5.0, random forests (RF), and elastic net (glmNet) classifiers were trained as presented in Figure 9 for prediction of any lung CT abnormalities at the 180-day follow-up visit. Variable importance measures (C5.0: % attribute/variable usage in the tree model (A); RF: difference in Gini index (B); glmNet: absolute value of the regression coefficient β (C)) for the 10 most influential explanatory variables are presented. CKD: chronic kidney disease; Q1, Q4: first, fourth quartile of anti-S1/S2 IgG titer; PD: pulmonary disease; CKD: chronic kidney disease.

Figure 9—figure supplement 5.. Variable importance statistics for prediction of moderate-to-severe lung computed tomography (CT) abnormalities at the 180-day follow-up by machine learning classifiers.

C5.0, random forests (RF), and elastic net (glmNet) classifiers were trained as presented in Figure 9 for prediction of moderate-to-severe lung CT abnormalities (severity score > 5) at the 180-day follow-up visit. Variable importance measures (C5.0: % attribute/variable usage in the tree model (A); RF: difference in Gini index (B); glmNet: absolute value of the regression coefficient β (C)) for the 10 most influential explanatory variables are presented. PD: pulmonary disease; GITD: gastrointestinal disease; Q1, Q2, Q4: first, second, fourth quartile of anti-S1/S2 IgG titer.

Figure 9—figure supplement 6.. Variable importance statistics for prediction of functional lung impairment at the 180-day follow-up by machine learning classifiers.

C5.0, random forests (RF), and elastic net (glmNet) classifiers were trained as presented in Figure 9 for prediction of functional lung impairment at the 180-day follow-up visit. Variable importance measures (C5.0: % attribute/variable usage in the tree model (A); RF: difference in Gini index (B); glmNet: absolute value of the regression coefficient β (C)) for the 10 most influential explanatory variables are presented. CKD: chronic kidney disease; Q1, Q2: first. second quartile of anti-S1/S2 IgG titer.

Figure 9—figure supplement 7.. Variable importance statistics for prediction of persistent symptoms at the 180-day follow-up by machine learning classifiers.

C5.0, random forests (RF), and elastic net (glmNet) classifiers were trained as presented in Figure 9 for prediction of persistent symptoms at the 180-day follow-up visit. Variable importance measures (C5.0: % attribute/variable usage in the tree model (A); RF: difference in Gini index (B); glmNet: absolute value of the regression coefficient β (C)) for the 10 most influential explanatory variables are presented. CVD: cardiovascular disease; GITD: gastrointestinal disease; COPD: chronic obstructive lung disease.

Figure 10.. Performance of the machine learning ensemble classifier in mild-to-moderate and severe-to-critical COVID-19 convalescents.

The machine learning classifier ensemble (Ens) was developed as presented in Figure 9. Its performance at predicting outcome variables at the 180-day follow-up visit (any computed tomography [CT] lung abnormalities, moderate-to-severe lung CT abnormalities [severity score > 5], functional lung impairment, and persistent symptoms) in the entire cohort, mild-to-moderate (outpatient or hospitalized without oxygen), and severe-to-critical COVID-19 convalescents (oxygen therapy or ICU) in repeated 20-fold cross-validation (five repeats) was assessed by receiver-operating characteristic (ROC) (Appendix 1—table 6). ROC curves and statistics (AUC: area under the curve; Se: sensitivity; Sp: specificity) in the cross-validation are shown. Numbers of complete observations and outcome events are indicated in the plots.