Radioproteomics stratifies molecular response to antifibrotic treatment in pulmonary fibrosis

Antifibrotic therapy with nintedanib is the clinical mainstay in the treatment of progressive fibrosing interstitial lung disease (ILD). High-dimensional medical image analysis, known as radiomics, provides quantitative insights into organ-scale pathophysiology, generating digital disease fingerprints. Here, we performed an integrative analysis of radiomic and proteomic profiles (radioproteomics) to assess whether changes in radiomic signatures can stratify the degree of antifibrotic response to nintedanib in (experimental) fibrosing ILD. Unsupervised clustering of delta radiomic profiles revealed 2 distinct imaging phenotypes in mice treated with nintedanib, contrary to conventional densitometry readouts, which showed a more uniform response. Integrative analysis of delta radiomics and proteomics demonstrated that these phenotypes reflected different treatment response states, as further evidenced on transcriptional and cellular levels. Importantly, radioproteomics signatures paralleled disease- and drug-related biological pathway activity with high specificity, including extracellular matrix (ECM) remodeling, cell cycle activity, wound healing, and metabolic activity. Evaluation of the preclinical molecular response–defining features, particularly those linked to ECM remodeling, in a cohort of nintedanib-treated fibrosing patients with ILD, accurately stratified patients based on their extent of lung function decline. In conclusion, delta radiomics has great potential to serve as a noninvasive and readily accessible surrogate of molecular response phenotypes in fibrosing ILD. This could pave the way for personalized treatment strategies and improved patient outcomes.


Supplementary Figures
Supplementary Figure 1.(A) Evaluation of the radiomic feature stability against inter-and intra-reader variation in the semi-automated lung segmentation workflow.Displayed is a representative transversal microCT image of a bleomycin-induced mouse lung.Outlined are three semi-automatically delineated lung contours of two different examiners (examiner 1: red and blue; examiner 2: yellow).For intra-and inter-operator intraclass correlation coefficient (ICC) analysis, a total of n=16 randomly selected lung scans covering different time points were segmented in this manner.(B) Boxplots displaying the distribution of the ICC coefficient per radiomic feature class for inter-operator ICC analysis and (C) intraoperator ICC analysis.The red dashed line indicates the set ICC threshold at 0.75.The stacked bar charts summarize the relative frequency and total number of robust and non-robust radiomic features.Supplementary Figure 2. Quantitative PCR of pro-fibrotic (Col1a1, Col3a1, Fn1), pro-inflammatory (Il6, Ccl2, Spp1), and nintedanib-targeted (Tgfb1, Timp1, Cxcl1, Ifng, Il1b, Tnf, Cd40lg) genes.Displayed is the mRNA fold change expression using the delta-delta Ct method (2 -ΔΔCt ) in cluster 1 (blue) and cluster 2 (red) compared to vehicle-treated samples.Each data point represents the mean of two technical replicates.Unpaired Student's t-test was used to compare the groups.Each entry is described by delta radiomic feature name, p-value, fold change difference, signed log10 enrichment p-value, cell type, and enrichment trend.Data provided in supplemental tables file.  1.

MicroCT image acquisition
Lung microCT images of each mouse were acquired at days 7 and 21 on a SkyScan 1176 (Bruker, Kontich, Belgium) in free-breathing conditions under isoflurane anesthesia using respiratory gated image acquisition.Anesthesia was induced by 5.0% and maintained by 1.5-2.5% isoflurane in air at 0.8-1.0L/min flow rate to achieve a breathing rate of 0.7-0.9breaths/s for an average scan time of 15 min.Animals were placed in supine position on the scanner bed with a styrofoam block mounted on the diaphragm to allow monitoring of respiratory gating.Image acquisition was performed with the following acquisition settings: tube voltage = 50 kV, tube current = 500 µA, filter = Al 0.5 mm, frame averaging = on (3), rotation step = 0.7 degrees, sync with events = 50 ms, X-ray tube rotation = 360 degrees, exposure time = 77 ms, resolution = 35 µm, slice thickness = 35 µm.Images were reconstructed with NRecon software (v.1.6.8.0; Bruker) using Feldkamp filtered back-projection algorithm with the following parameters: misalignment compensation (scan-dependent manual adjustment), smoothing = 1 with Gaussian kernel, ring artifact compensation = 4, and beam hardening correction = 10%.Reconstructed images were converted to DICOM format.

CT segmentation of mouse lungs
Left and right lung lobes of mice were semi-automatically segmented by two readers (D.L., M.B.) using MIM software (v.7.1.6,MIM Software Inc., Cleveland, Ohio, USA).Briefly, a seed was set within the right and left lung using the "region grow" tool (upper limit = -600 HU, lower limit = -800 HU, tendril diameter = 0.2 mm, fill holes = strong), which then automatically defined the vast majority of the lung volume in the 3D space.Manual contour alignment with the 2D/3D brush was used to correct misaligned areas.Finally, the smoothing function was used to remove sharp edges.For medical diagnostics, the Hounsfield scale is usually normalized to 120 kVp tube voltage, which could technically not be achieved by our microCT instrument.To enable direct comparison between CT-derived radiomic datasets from patients and mice, the reconstructed microCT images of mouse lungs were pixel value corrected to match clinical specifications as previously described (10).
HRCT image acquisition HRCT acquisition of human lungs was performed at Bern University Hospital or outpatient clinics.Instrument and scan settings used are summarized in Supplementary Table 11.All HRCT scans were evaluated by a senior radiologist (L.E.) at the Department of Diagnostic, Interventional, and Pediatric Radiology of the Bern University Hospital for the presence of PF-ILD on a standard picture archiving and communication system workstation and a radiology-grade display monitor.

CT segmentation of human lungs
Left and right lung lobes were semi-automatically segmented by two readers (C.M., L.K.) with the opensource software 3D Slicer (v.5.2.1).Pulmonary hilar vessels and atelectatic areas were manually excluded from the regions of interest.Manual contour corrections were only applied when spatially limited areas did not coincide with the actual borders of the lungs.

Pulmonary function tests
Pulmonary function tests were performed by trained personnel at the Department of Pulmonary Medicine of the Bern University Hospital or in outpatient clinics.All tests were performed following established protocols (11)(12)(13)(14).
Radiomic feature calculation Calculation of radiomic features was performed on merged structures of left and right lung lobes using Z-Rad software (v.7.3.1, https://medical-physics-usz.github.io/,Department of Medical Physics, University Hospital Zurich, Zurich, Switzerland), an image biomarker standardization initiative (IBSI)compliant Python-based software (15), as described in (10).Mouse lungs were resized to isotropic voxels of 0.15 mm.To achieve comparable voxel size in patients, human lungs were resized to isotropic voxels of 2.75 mm, corresponding to an estimated 6000-fold volumetric difference (16).Both mouse and human lung volumes were discretized to a fixed bin size of 50 HU in a range of -1000 HU to 200 HU.From the resized volumes, 1388 radiomic features were calculated per lung scan and time point (HU limits: -1000 to 200 HU), corresponding to the following feature classes: 1 Histogram features carry information about distribution of voxel intensities using first-order statistics (e.g.mean, standard deviation, skewness, kurtosis), describing tissue intensity characteristics.Texture features define intra-tissue heterogeneity by calculating the spatial relationship between neighboring voxel intensities (17).Wavelet features compute histogram and texture features after wavelet decompositions of the original image using eight different coiflet filters (high-to low-pass filters), thereby concentrating the features on different frequency ranges (18).Shape features describe tissue volume and size independent of intensity distribution.Delta radiomic features describing the change of each feature between pre-and post-treatment were expressed as delta values: ΔFeature = Feature (t2) -Feature (t1) (19).Lung densitometric information was directly inferred from the radiomic histogram feature hist_mean, which describes the lung attenuation-based average HU intensity of the segmented lung volume.
Radiomic feature stability evaluation Intraclass correlation coefficients (ICC) were calculated for each radiomic feature to evaluate stability against inter-and intra-reader bias in the lung segmentation process.For inter-and intra-reader ICC, two (D.L., M.B.) and one examiner(s) (D.L.), respectively, independently segmented 16 randomly selected mouse lung scans, followed by radiomic feature calculation of the delineation structures.ICCs were calculated using two-way mixed effect models with the consistency method in the irr R package according to published reports (20).Only stable/reproducible features (n=1130) with ICC≥0.75 were considered for further analyses for both mouse and human datasets (21).Feature stability assessment was performed on the mouse dataset due to the lesser degree of automation in lung segmentation.

Proteomics
For comparative proteomics, the middle lobe of the right mouse lung was snap frozen in liquid nitrogen and stored at -80°C until processing.Sample workup and data collection was performed by trained personnel at the Proteomics and Mass Spectrometry Core Facility (PMSCF) at the University of Bern using standard established protocols.All vehicle-(n=14) and nintedanib-treated (n=10) samples were analyzed.One vehicle sample was excluded from analysis due sample workup issues.Tissue homogenization was performed in 8M urea / 100 mM Tris (pH 8.0) buffer supplemented with cOmplete protease inhibitor cocktail (Roche Diagnostics, Mannheim, Germany) using the FastPrep system (MP Biomedicals).Following reduction, alkylation, and overnight protein precipitation with ice-cold acetone, 10 µg of the cleaned protein mixture was digested into peptides using a two-step digestion protocol (LysC for 2 h at 37°C followed by Trypsin at room temperature overnight).Digests were analyzed by nano-liquid chromatography on a Dionex Ultimate 3000 (ThermoFisher Scientific, Reinach, Switzerland) through a CaptiveSpray source (Bruker, Bremen, Germany) with an end-plate offset of 500 V, a drying temperature of 200°C, and with the capillary voltage fixed at 1.6 kV.A volume of 2 µL (200 ng) protein digest was loaded onto a pre-column (PepMap 100 C18, 5 µm, 100 A, 300 µm diameter x 5 mm length, ThermoFisher) at a flow rate of 10 µL/min with 0.05% trifluoroacetic acid in water / acetonitrile 98:2.After loading, peptides were eluted in back flush mode onto an in-house made C18 CSH Waters column (1.7 µm, 130 Å, 75 µm x 20 cm) by applying a 90-minute gradient of 5% acetonitrile to 40% in water / 0.1% formic acid, at a flow rate of 250 nL/min.The timsTOF Pro instrument (Bruker, Bremen, Germany) was operated either in data-dependent acquisition (DDA) or data-independent (DIA) mode using the Parallel Acquisition Serial Fragmentation (PASEF) option.The mass range was set between 100 and 1700 m/z, with 10 PASEF scans between 0.7 and 1.4 V s/cm 2 .The accumulation time was set to 2 ms, and the ramp time was set to 100 ms, respectively.Fragmentation was triggered at 20'000 arbitrary units, and peptides (up to charge of 5) were fragmented using collision induced dissociation with a spread between 20 and 59 eV.DDA data was processed further with FragPipe software (v.17.0) using the IonQuant algorithm and filtering protein identifications to a 1% false discovery rate (FDR) on the peptide level using the Percolator algorithm.Furthermore, protein groups were filtered by the criterion that at least two different razor peptide sequences were identified as evidence for the existence of the protein group.From the DDA data, a spectral library was built with the FragPipe software.This library was used to identify and quantify proteins with the DIA data using standard parameters in Spectronaut 16 software (Biognosys, Schlieren, Switzerland).Protein names (Uniprot IDs) were converted to Entrez IDs and Gene Symbols using UniProt.ws and annotationDbi R packages.Protein names without matching Entrez Gene ID were dropped, resulting in a final set of 7006 proteins.

Phosphoproteomics
For phosphoproteomics, the middle lobe of the right mouse lung was snap frozen in liquid nitrogen after collection and stored at -80°C until processing.Sample workup and data pre-processing was performed by trained personnel at the Proteomics and Mass Spectrometry Core Facility (PMSCF) at the University of Bern using standard established protocols.Randomly selected subsets of vehicle-(n=5) and nintedanib-treated (n=5) were analyzed.A titanium dioxide phosphopeptide enrichment workflow (22) with subsequent DDA liquid chromatography tandem mass spectrometry (LC-MS) analysis on the same instrument and parameter settings as described above was applied.Samples were searched and quantified with FragPipe (23) (v.18.0,MSFragger version 3.5, Philosopher version 4.4.0,IonQuant version 1.8.0) using the following parameters: swissprot (24) Mus musculus database (release 2022_01) with isoforms and common contaminants; 20 ppm and 0.05 Da mass tolerance for precursors and fragment, respectively; search enzyme trypsin with max 3 allowed missed cleavages; fix modification: carbamidomethylation of cysteine; variable modifications (altogether max 4/peptide): methionine oxidation of methionine (max 3/peptide), phosphorylation of serine, threonine and tyrosine (max 3/peptide) and protein N-terminal acetylation.Peptide forms normalized with the variance stabilization (25) normalization method are reported as NormI, along FragPipe's MaxLFQ and IonQuant's FragI abundance measures.The intensities of peptide forms were combined as protein phosphosite locations by summing the corresponding contributions.
Kinase activity enrichment analysis Differential expression of phosphosites and subsequent kinase activity enrichment analysis was performed as described in (22).First, the phosphosites' missing values were imputed using a leftcensored Gaussian replacement method if there was more than 1 missing value in a group of replicate, and a maximum likelihood estimation otherwise (26).A moderated t-statistic (26) was then calculated for each phosphosite, and used as the ranking metric for the Kinase Activity Enrichment Analysis (KAEA) tool (22).KAEA was then applied on the ranked phosphosite list and reversed ranked list using the included mouse kinase substrate database.SetRank set p-value and FDR cutoff were set to 0.01 and 0.05, respectively.Differential protein expression analysis Differential expression of proteins between groups of interest was calculated in R using the "limma" package according to standard guidelines (27).At first, DIA-based Spectronaut protein expression intensities were log2-transformed.Then, log2 fold changes were calculated as contrasts by application of a linear model using robust regression for each protein.Finally, estimated coefficients and standard errors for the given set of contrasts were calculated for each protein, followed by Empirical Bayes smoothing of standard errors.Proteins with log2FC>0.3(p<0.05),corresponding to 23% mean expression change, were considered as statistically significant.
Gene expression analysis RNA was isolated from blood-free cranial lobes of the right mouse lung stored in RNAlater (ThermoFisher Scientific).Tissues were mechanically homogenized with the TissueLyser II instrument (Qiagen, Hombrechtikon, Switzerland), followed by total RNA isolation with the RNeasy Tissue Mini Kit (Qiagen, Hombrechtikon, Switzerland).Isolated RNA was reverse transcribed into cDNA using the Transcriptor First Strand cDNA Synthesis Kit (Roche Diagnostics, Switzerland).Expression of fibrotic (Col1a1, Col3a1, Fn1), inflammatory (Il6, Ccl2, Spp1), and nintedanib-related (Tgfb1, Timp1, Cxcl1, Ifng, Il1b, Tnf, Cd40lg) genes was analyzed by SYBR Green quantitative PCR using GoTaq Green Master Mix kit (Promega) as described in (28).Expression of mRNA was expressed to delta Ct values (Ct [gene of interest] -Ct [reference gene]) with Rplp0 as reference gene.Lower delta Ct values indicate higher target gene expression.Fold changes relative to vehicle-treated samples were calculated using the delta-delta Ct method.The list of the primer pairs used in this study is provided in Supplementary Table 12.

Immunofluorescence and microscopy
Formalin-fixed paraffin-embedded lung sections (3 µm thickness) were cut on a HistoCore Multicut microtome (Biosystems Switzerland AG, Muttenz, Switzerland).Following deparaffinization, heatmediated antigen retrieval with R-Universal Buffer (Cat.AP0530-500, Aptum Biologics) was performed for 15 min at 95°C.After incubation for 25 min at RT for cooling, blocking of unspecific antibody staining was performed with 5% BSA in antibody diluent (Cat.S3022, Dako) for 1 h at RT.Primary antibodies dissolved in antibody diluent (Cat.S3022, Dako) were then applied and incubated overnight at 4°C.Next, samples were incubated with secondary antibodies dissolved in PBS supplemented with 1% BSA for 2 h at RT.All antibodies and the dilutions used are listed in Supplementary Table 13.Finally, cell nuclei were counterstained with 4',6-diamidino-2-phenylindole (DAPI) for 10 min at RT.The sections were then scanned in immunofluorescence mode on a AxioScan.Z1 slide scanner (Zeiss, Feldbach, Switzerland) using a Plan-Apochromat 20x/0.8M27 objective.Cells positively stained for α-SMA were quantified using the "Positive cell detection" tool of the open source software QuPath (v.0.4.0) at default parameter settings and detection thresholds 800, respectively (29).From each sample, five representative areas at 500x500 µm were quantified and the average was used for statistical analyses.

Unsupervised clustering
All variables were z-scored ([x-mean]/standard deviation) prior to analysis (clusterSim R package).Unsupervised agglomerative hierarchical (using Euclidean distance with complete linkage method) or k-means clustering was performed to identify subgroups of mice or patients with similar delta radiomic feature patterns or proteomic profiles using base R functions.Clusterability was evaluated by Hopkin's statistic H, with H>0.5 indicating clusterability (hopkins R package) (30).The optimal number of clusters was determined by average silhouette statistics inspecting k clusters between 2 and 5, selecting the optimal k based on global or local optimum for separation (factoextra R package).Stability of clusters was assessed by Jaccard bootstrapping with n=1'000 iterations (fpc R package) (31).

Variable importance evaluation
To estimate the importance of each delta radiomic feature for the classification produced by unsupervised clustering, we calculated filter-based variable importance using the "caret" package (32) and retained features with classification score≥0.9.Features (n=54) most important for differentiating clusters 1 and 2 in nintedanib-treated mice are listed in Supplementary Table 14.
Gene Ontology and Reactome pathway enrichment Curated lists of DE or highly-correlated proteins were used to perform Gene Ontology (GO) or Reactome pathway enrichment analysis using the R packages clusterProfiler and ReactomePA, respectively (33)(34)(35), retaining results after false discovery rate adjustment (p<0.05).In case of enrichment of proteins highly correlated with delta radiomic features, proteins with positive and negative correlation coefficients were entered separately into GO or Reactome pathway enrichment analysis.To visualize and interpret results, the GeneRatios of positive and negative enriched pathways were transformed into matrices with delta radiomic features as columns and pathways terms as rows.Then, results with GeneRatios<0.10were dropped (set to zero), and only pathway terms enriched (GeneRatio≥0.10) in at least two delta radiomic features were retained.Subsequently, the tables containing positive and negative enriched delta radiomic pathway pairs were aggregated, and rows and columns without significantly enriched results were removed, followed by visualization with the pheatmap R package.
Correlation analysis Spearman's rank correlation coefficient rho was calculated between selected delta radiomic features and the log2-transformed expression intensity of every protein using base R packages, retaining only proteins with p<0.05 and ⍴≥0.6 for further analysis.Pearson's correlation coefficient r was calculated between selected delta radiomic features and the fraction of α-SMA positive cells using base R packages.
Cell type signature enrichment analysis To infer relative cell type frequency changes between two groups from proteomics data, we applied signature enrichment analysis as described in (36,37), utilizing their published dataset.Cell type signatures were defined as sets of genes with cell-type specific gene expression of log2 fold change>0.3and adjusted p<0.05 (37).For each cell type, we then tested for the enrichment in a ranked list of DE proteins (log2 fold changes) or correlation coefficients (weighted by -log10 p-value) using the Kolmogorov-Smirnov test.Positive and negative signed enrichment scores (-log10 p-values signed by effect size) reflect relative depletion and enrichment of the respective cell type, respectively.To visualize and interpret cell type signatures enriched in the sets of proteins highly-correlation (Spearman's |⍴|≥0.6,p<0.05) with delta radiomic features, the signed enrichment scores of each variable were transformed into a matrix with delta radiomic features as columns and cell types as rows.Enriched results with enrichment score<2 (corresponding to p<0.01) were dropped (set to zero), followed by removal of rows and columns without enriched entries, and visualization with the pheatmap R package.
Association analysis with clinical parameters Association analyses were performed to investigate associations of patient delta radiomics-derived (kmeans) clusters with clinical parameters.Mann-Whitney U test was used for comparison of numerical variables, and Fisher's exact test was used to compare categorical variables.

Statistical analyses
All statistical analyses were performed in R (v. 4 K-means clustering of z-scored delta radiomic features (n=244) of nintedanib-treated mice (n=10) shows two stable clusters (Jaccard coefficients>0.90,where 1 describes perfect stability).(E) Heatmap summary of the k-means clustering results (nintedanib-treated mice, n=10).Clusters and the feature class of each variable are indicated.(F) Lung tissue density in cluster 1 and 2 expressed as mean Hounsfield unit (HU) intensity post-treatment.Mann-Whitney U-test was used to compare groups.(G) Principal component analysis of the phosphosite expression values (n=20043) in subsets of randomly selected nintedanib-(n=5) and vehicle-treated (n=5) mice.(H) Kinase activity enrichment analysis (KAEA) of differentially expressed phosphosites in nintedanib-against vehicle-treated mice.Underactive kinases are colored in red, over-active kinases are colored in blue.(I) Volcano plots of protein expression in cluster 1 and (J) cluster 2 compared to vehicle-treated mice.Proteins with log2FC>0.30and p<0.05 were considered to be differentially expressed.Down-and upregulated proteins are highlighted in blue and red, respectively.

Supplementary Figure 4 .
each variable are indicated.(D) Heatmap displaying Gene Ontology Biological Process (GO:BP) pathways enriched (GeneRatio≥0.10,adjusted p<0.05) in radioproteomic association modules for positively (Spearman's ≥0.6,p<0.05, red annotation) or negatively (Spearman's ≤-0.6,p<0.05, blue annotation) correlating proteins.Only pathways enriched in at least two radioproteomic association modules are displayed.Association modules without enriched pathways following filtering are not displayed.(E) Heatmaps displaying the results of unsupervised k-means clustering of z-scored subsets of delta radiomic features positively enriched in extracellular matrix (ECM) remodeling (n=8) and wound healing (n=7) in nintedanib-treated mice (n=10), respectively.Cluster assignment of samples and Reactome pathway enrichment of variables are indicated.Delta radiomics stratifies the degree of pulmonary function decline in nintedanib-treated PF-ILD patients.(A) Left: K-means cluster plot for z-scored preclinical responsedefining delta radiomic features (n=54) of the PF-ILD cohort (n=19) indicates three fairly robust clusters (Jaccard coefficients>0.60,where 1 describes perfect stability).Right: Scatter plot showing the average silhouette coefficient versus the number of clusters for the k-means clustering input data.The blue dashed line indicates the global optimum.(B) Box plots comparing age (years) # and disease duration (years) ¶ at baseline between clusters A1-A3.Mann-Whitney U test was used to compare the groups.(C) Associations of clusters A1-A3 with clinical and demographic parameters in the PF-ILD cohort.Fisher's exact test was used to compare the categorical variables.(D) Left: K-means cluster plot for zscored delta radiomic features positively correlating with ECM-remodeling (n=8) of the PF-ILD cohort (n=19) indicates two stable clusters (Jaccard coefficients>0.75,where 1 describes perfect stability).Right: Scatter plot showing the average silhouette coefficient versus the number of clusters for the kmeans clustering input data.The blue and red dashed line indicate the global and local optimum, respectively.(E) Box plots comparing age (years) # and disease duration (years) ¶ at baseline between clusters B1 and B2.Mann-Whitney U test was used to compare the groups.# Age is defined as the period between birth date and baseline HRCT scan.¶ Disease duration is defined as the period between the first reported diagnosis of PF-ILD and baseline HRCT scan.(F) Associations of clusters B1 and B2 with clinical and demographic parameters in the PF-ILD cohort.Fisher's exact test was used to compare the categorical variables.

Table 9 .
(1)ociations of patients' groups resulting from unsupervised clustering on preclinical treatment response-defining delta radiomic features (n=54) with clinical parameters.ILD, interstitial lung disease; SSc, systemic sclerosis; HP, hypersensitivity pneumonitis; FVC, forced vital capacity; FEV1, forced expiratory volume in 1 s; DLCO, diffusing capacity of the lung for carbon monoxide; P(A)H, pulmonary (arterial) hypertension.NA denotes missing values as n (%).*:Age at time of baseline (pre-treatment) HRCT scan.†:diseaseduration was defined as the period (months) between first reported diagnosis of PF-ILD in the patient records and the baseline (pre-treatment) HRCT scan.‡:PH was assessed by echocardiography or right heart catheterization.If right heart catheterization was performed, mPAP>20 mmHg was considered diagnostic(1).§: Concomitant immunomodulatory therapy was indicated if a patient received immunomodulatory medication at any time simultaneous to nintedanib treatment.Immunomodulatory therapy included prednisolone, mycophenolate mofetil, azathioprine, rituximab, tocilizumab, or combinations thereof.
Data are presented as median (± interquartile range (IQR)) or n (%).Mann-Whitney U and Fisher's exact tests were used to compare the numerical and categorical variables, respectively.Abbreviations: IPF, idiopathic pulmonary fibrosis;

Table 10 .
(1)ociations of patients' groups resulting from unsupervised clustering on ECM remodeling-associated delta radiomic features (n=8) with clinical parameters.Data are presented as median (± interquartile range (IQR)) or n (%).Mann-Whitney U and Fisher's exact tests were used to compare the numerical and categorical variables, respectively.Abbreviations: IPF, idiopathic pulmonary fibrosis; ILD, interstitial lung disease; SSc, systemic sclerosis; HP, hypersensitivity pneumonitis; FVC, forced vital capacity; FEV1, forced expiratory volume in 1 s; DLCO, diffusing capacity of the lung for carbon monoxide; P(A)H, pulmonary (arterial) hypertension.NA denotes missing values as n (%).*:Age at time of baseline (pre-treatment) HRCT scan.†:diseaseduration was defined as the period (months) between first reported diagnosis of PF-ILD in the patient records and the baseline (pre-treatment) HRCT scan.‡:PH was assessed by echocardiography or right heart catheterization.If right heart catheterization was performed, mPAP>20 mmHg was considered diagnostic(1).§: Concomitant immunomodulatory therapy was indicated if a patient received immunomodulatory medication at any time simultaneous to nintedanib treatment.Immunomodulatory therapy included prednisolone, mycophenolate mofetil, azathioprine, rituximab, tocilizumab, or combinations thereof.

Table 11 .
Summary of HRCT scan acquisition parameters.

Table 12 .
Mouse primer sequences used for quantitative PCR reactions.

Table 13 .
Primary and secondary antibodies used for immunofluorescence staining.

Table 14 .
Variable importance selected delta radiomic features used to classify treatment response. .Pre-treatment PFT (max. 1 month after treatment initiation) b.Post-treatment PFT (min.6 months interval to pre-treatment PFT) Exclusion Criteria 5. Presence of secondary lung disease at times any HRCT scans and PFT recording (e.g., cancer, COVID-19, pneumonia, bronchitis) Demographic and clinical parameters were derived from electronic patient records, including age (birth date), sex, disease etiology, date of diagnosis, date of nintedanib treatment start, presence of pulmonary (arterial) hypertension, smoking status, concomitant medications, and dates of PFT and HRCT scan recordings.The recorded PFT parameters included forced vital capacity (FVC) in % pred.and liters.Changes in PFT recordings between pre-and post-treatment were expressed as delta values.A Summary of patient demographics and clinical characteristics is provided in Table Abbreviations: GLCM = Gray Level Co-occurrence Matrix, NGTDM = Neighborhood Gray Tone Difference Matrix, GLRLM = Gray Level Run Length Matrix, GLDZM = Gray Level Distance Matrix and NGLDM = Neighboring Gray Level Dependence Matrix.a .3.1.)environment.For all analyses, a p<0.05 was considered statistically significant unless stated otherwise.The following R packages were used: AnnotationDbi, caret, clusterProfiler, clusterSim, doParallel, dplyr, factoextra, fpc, ggplot2, ggpubr, ggrepel, limma, openxlsx, org.Mm.eg.db, parallel, pheatmap, readxl, rstatix, tidyverse, UniProt.ws,VennDiagram.Visualization Figures were created in Adobe Illustrator (v.28.2) and partially contain graphics or illustrations from Adobe Stock and BioRender.comaccessed through the academic licenses of the University of Bern.