The study population consisted of hospitalized patients from two centers: A derivation cohort from April 2014 to October 2019 and a longitudinal internal validation cohort from November 2019 to May 2020 were recruited in Renji Hospital (Shanghai). However, an external validation cohort from Guangzhou Institute of Respiratory Health was constructed during the same time period. Ethical approval was obtained for this retrospective analysis, and the need to obtain informed consent was waived for both institutions. Eligible patients for this study initially fulfilled Bohan and Peter's criteria for DM or Sontheimer's criteria for clinically amyopathic DM on admission (14, 15) and were retrospectively confirmed by the 239th European NeuroMuscular Center (ENMC) classification criteria for DM (16). All patients presented with imaging-confirmed ILD and positive anti-MDA5 antibody. The anti-MDA5 antibody was detected by immunoblotting assay (Euroimmun, Germany) and confirmed by enzyme-linked immunosorbent assay (Supplementary Material). ILD course was defined as time from the first abnormal pulmonary CT, which revealed ILD changes to admission. Patients with ILD course > 3 months or coexisting malignancy (within 3 years) or preexisting chronic obstructive pulmonary disease were excluded. The primary outcome was the 6-month all-cause mortality since the time of admission.

The derivation cohort was randomly divided into the training (n = 121) and testing (n = 31) sets. The internal and external validation cohorts comprised 44 and 32 patients, respectively (Figure 1).

Patients underwent non-contrast pulmonary HRCT within a week of admission (median, 2 days; range, 1–6 days) using a multidetector CT scanner. CT slice thickness was 1.0–1.5 mm at 10-mm interval with standard reconstruction. Radiomics analysis was performed using a dedicated software (Radiomics, syngo.via Frontier 1.2.2, Siemens Healthineers, Forchheim, Germany). This platform integrated an open-source library named “PyRadiomics” (version 2.0.1). First, the Digital Imaging and Communications in Medicine files of CT images were imported; segmentation of the volumetric lung was performed automatically with deep learning algorithm and adjusted manually by ZC and YZ (both with more than 5 years' experience in chest CT imaging). The medical image series were resampled to 1 1 1 mm voxel size with B-spline interpolator before subsequent feature extraction steps. Bin width was set at 25 to create histograms for discretization of the image gray levels. The extracted features are as follows: (i) 18 first-order features, which describe the distribution of voxel intensities in region of interest (ROI); (ii) 18 shape and size features, which include descriptors of 2D and 3D size and shape of ROI; (iii) 75 texture features, describing patterns or the spatial distribution of voxel intensities. In addition to the original images, wavelet-transformed images were generated using eight different high- to low-frequency band combinations according to three directions.

Thus, a total of 855 features, namely, 111 from the original image and 744 from the wavelet-transformed images, were outputted. Then, the intra- and inter-observer reliabilities were evaluated using Spearman's correlation coefficients in derivation cohort; features with rho>0.8 were included for further analysis.

The radiomics model was derived from the training set and tuned in the testing set. The least absolute shrinkage and selection operator (LASSO) and Cox regression model was used to select features, with penalty parameter tuning conducted by ten-fold cross-validation (17). Then, optimal features for radiomics model were determined by backward stepdown analysis with Akaike's information criterion (AIC) applied as the stopping rule.

Radiomic score for each patient was calculated from a linear combination of selected radiomic features weighted by respective coefficients.

All CT images were reviewed by ZC and YZ who were blinded to patients' outcome. Inter-observer variability was evaluated by intraclass correlation coefficient (ICC). The results were agreed upon by consensus between the two observers. HRCT findings were graded on a scale of 1–6 based on the classification system that corresponded to consecutive pathologic phases as previously reported (8): 1, normal attenuation; 2, ground-glass attenuation (GGA); 3, consolidation; 4, GGA associated with traction bronchiolectasis or bronchiectasis; 5, consolidation associated with traction bronchiolectasis or bronchiectasis; and 6, honeycombing.

The lungs were divided into six zones (upper, middle, and lower on both sides); each zone was evaluated separately. The score was based on the percentage of the abnormal lung parenchyma and was estimated to the nearest 5% of parenchymal involvement. The overall CT visual score was calculated by adding the average score of six zones (6).

In addition, an ILD-GAP score without lung imaging information was also calculated in our analysis as a comparator, which was modified from the IPF-derived GAP staging system and had been validated in some non-IPF ILD cohorts for prognosis prediction (18, 19).

Finally, a subgroup analysis was performed to show in which cases the radiomics model outperformed the visual scoring model.

The differences of clinical parameters between the derivation and internal validation and external validation cohorts were assessed by the Mann–Whitney U-test, chi-square test, or Fisher exact test, as appropriate. Post hoc comparisons were made with Bonferroni's correction for adjustments of multiple comparisons. The optimal cutoff values for significant continuous variables were identified by X-tile (20). The association between Rad-score and outcome was assessed by Kaplan–Meier survival plot and log-rank test.

An improved multidimensional prediction model based upon Rad-score (Rad-score plus model) was constructed by integrating clinically significant predictors from univariable analysis into the multivariable Cox proportional hazards model (backward stepdown selection with AIC).

Correlation analysis between visual score, Rad-score, and other variables were applied with Spearman's coefficients. Model discrimination was quantified by the Harrell concordance index (C-index) with 95%CI. The change in overall log-likelihood ratio was used to assess the increase in predictive power. Additionally, the model calibration was assessed using the calibration curve, Greenwood–Nam–D'Agostino (GND) test, and Brier score (13, 21, 22). A decision curve analysis was built to determine and compare the clinical usefulness of each model (23). Significance was defined as p-value <0.05.

Statistical analyses were performed on SPSS software version 25 (IBM Corp., Armonk, NY, United States) and R software version 3.6.1. All R codes used were available at GitHub (https://github.com/tomato08217/MDA-ILD5).

Eight hundred and twenty-five features remained after reproducibility analysis. Five radiomic features were finally selected for modeling as follows: Original_Shape_Flatness (A), Wavelet.LLL_Glcm_MCC (B), Wavelet.LLH_Glcm_Ldmn (C), Wavelet.HLL_Glszm_LargeAreaHighGrayLevelEmphasis (D), and Wavelet.LLL_Firstorder_Skewness (E).

The contribution of the selected parameters with their regression coefficients was presented in the form of a histogram in Figure 1. Then, the Rad-score was computed using the following formula: A*(0.599) +B*(0.461) +C*(0.408) +D*(−0.543) +E*(−1.19). The C-index values of the radiomics model were 0.83 (95%CI, 0.74–0.91), 0.82 (95%CI, 0.65–0.99), 0.78 (95%CI, 0.64–0.93), and 0.76 (95%CI, 0.59–0.92) for the training, testing, internal validation, and external validation sets, respectively. The GND test showed good calibration abilities in the four datasets. The prediction errors were low, attested by Brier score ranging from 0.12 to 0.19.

The optimal cutoff value for Rad-score was determined as 1, dividing patients into the high-risk group (Rad-score > 1) and low-risk group (Rad-score ≤ 1) with distinct 6-month mortality (Figure 2A). The representative patients' CT images of the two groups are presented in Figures 2B,C, respectively, with corresponding Rad-scores.

The univariable Cox regression results of baseline clinical characteristics and treatment between the survivors and deceased in the derivation cohort are presented in Table 2. Then, 10 candidate predictors, i.e., with significant differences in univariable analyses, were incorporated into a multivariable Cox proportional hazards model with Rad-score, indicated as follows: age on admission, course of DM, arthralgia, three-category FVC% (FVC% ≥ 50%, FVC% <50%, and unable to perform pulmonary function tests), arterial oxygen/fraction of inspiration oxygen (PaO₂/FiO₂), lactate dehydrogenase (LDH), serum ferritin (ng/mL), C-reactive protein, lymphocyte, and maximum dosage of methylprednisolone (Table 3).

The optimal combination obtained by the backward stepwise algorithm with minimum AIC was as follows: Rad-score, FVC% (three categories), and age (on admission). The respective hazard ratios and β-regression coefficients are presented in Table 4. The final Rad-score plus model was developed to predict the 6-month mortality with C-index values of 0.88 (95%CI, 0.79–0.96) in the training set, 0.88 (95%CI, 0.71–1.0) in the testing set, 0.83 (95%CI, 0.68–0.98) in the internal validation cohort, and 0.84 (95%CI, 0.64–1.0) in the external validation cohort, respectively. Nomogram was presented online to provide individualized risk estimates (https://mda5-ild.shinyapps.io/MDA_ILD5_test/).

The agreement of visual CT score between two readers was good with ICC of 0.821 (95%CI, 0.711–0.898). The visual score model yielded the C-index values of 0.75 (95%CI, 0.66–0.84) for the training set, 0.81 (95%CI, 0.65–0.98) for the testing set, and 0.79 (95%CI, 0.65–0.94) for the internal validation cohort, when predicting the 6-month mortality of MDA5⁺ DM-ILD. However, the discriminative power of visual score dramatically declined in the external validation cohort with the C-index of 0.66 (95%CI, 0.49–0.82) and poor calibration ability. However, the C-index values of ILD-GAP model were relatively even in each dataset (see Table 5 for details).

The prediction accuracy and robustness of the ILD-GAP score, visual score, radiomics model, and Rad-score plus model are compared in Table 5, with respective calibration curves (Figure 3). The final Rad-score plus model had the best performance surpassing the radiomics model, visual score, and ILD-GAP score (p < 0.001). Finally, the decision curve analysis demonstrated that, in terms of clinical applicability, the Rad-score plus model also had a higher overall net benefit than other models across the majority of the ranges of threshold probabilities (Figure 4).

In the correlation analysis, the visual score, Rad-score, and Rad-score plus mentioned above all showed strong correlation with FVC and PaO₂/FiO₂ rather than other clinical predictors (for details, please see Supplementary Table 3).

Across all datasets, radiomics performed better than visual score in 35 cases, which were defined as Rad-score better group. This subgroup presented a significantly higher mortality, a shorter disease course, a worse pulmonary function, and a higher dosage requirement of steroid (Supplementary Table 2), possibly implying that the radiomics model has a better performance in the patients with more progressive disease course.