Improved prediction of radiation pneumonitis by combining biological and radiobiological parameters using a data-driven Bayesian network analysis.

Grade 2 and higher radiation pneumonitis (RP2) is a potentially fatal toxicity that limits efficacy of radiation therapy (RT). We wished to identify a combined biomarker signature of circulating miRNAs and cytokines which, along with radiobiological and clinical parameters, may better predict a targetable RP2 pathway. In a prospective clinical trial of response-adapted RT for patients (n = 39) with locally advanced non-small cell lung cancer, we analyzed patients' plasma, collected pre- and during RT, for microRNAs (miRNAs) and cytokines using array and multiplex enzyme linked immunosorbent assay (ELISA), respectively. Interactions between candidate biomarkers, radiobiological, and clinical parameters were analyzed using data-driven Bayesian network (DD-BN) analysis. We identified alterations in specific miRNAs (miR-532, -99b and -495, let-7c, -451 and -139-3p) correlating with lung toxicity. High levels of soluble tumor necrosis factor alpha receptor 1 (sTNFR1) were detected in a majority of lung cancer patients. However, among RP patients, within 2 weeks of RT initiation, we noted a trend of temporary decline in sTNFR1 (a physiological scavenger of TNFα) and ADAM17 (a shedding protease that cleaves both membrane-bound TNFα and TNFR1) levels. Cytokine signature identified activation of inflammatory pathway. Using DD-BN we combined miRNA and cytokine data along with generalized equivalent uniform dose (gEUD) to identify pathways with better accuracy of predicting RP2 as compared to either miRNA or cytokines alone. This signature suggests that activation of the TNFα-NFκB inflammatory pathway plays a key role in RP which could be specifically ameliorated by etanercept rather than current therapy of non-specific leukotoxic corticosteroids.

Introduction

Chemoradiation (ChemoRT) followed by the immune checkpoint blocker (ICB) durvalumab has become the new standard of care for patients with locally advanced NSCLC patients [1]. Unfortunately, treatment of locally advanced NSCLC is complicated by the frequent occurrence of symptomatic (grade 2 or higher) radiation pneumonitis (RP) which may occur in up to 25% of patients [2], [3], [4]. Multiple recent reports noted a rise in symptomatic Grade 2+ pneumonitis with adjuvant ICB [5], [6], [7]. In a nationwide cohort of more than 2,000 patients, we found that grade 3+ pneumonitis was the most common toxicity in locally advanced lung patients [8] and that it led to discontinuation of potentially curative durvalumab in more than 10% of “real-world” patients as compared to 3% enrolled on PACIFIC trial (which only scored grade 3 RP and included only patients who successfully completed chemoRT without RP). High dose corticosteroids are given to symptomatic patients with Grade 2+RP, but this induces leukopenia which may negatively impact ICB efficacy [9,10]. We have found on multivariate analysis that the development of immune-related Adverse Events (irAE) in locally advanced NSCLC patients significantly reduced overall survival [8]. Therefore, there is an unmet clinical need to develop immune sparing RP treatment strategies.

Lung irradiation causes the release of inflammatory cytokines and miRNAs that can promote early acute RP and late pulmonary fibrosis (PF) [11,12]. Among various cytokines, early release of TNFα plays a priming role in RP [13]. Thus, targeting of TNFα signaling is a possible strategy to develop a lung radioprotector [14]. Our preclinical data showed that inhibition of TNFα signaling by genetically deleting the receptor 1 (TNFRSF1A-null mice) or via antisense oligonucleotide (ASO) protects the lung from radiation toxicity [14]. We further extended our observation using Tristetraprolin (TTP, Zfp36 in mice) knockout mice, which have high basal levels of TNFα. TTP is a key negative regulator of TNFα by causing its mRNA degradation. These mice develop RP symptoms within a week post-irradiation [15]. In the TNFα-TNFR1 mediated inflammatory signaling, TNFα converting enzyme (TACE), also known as A disintegrin and metalloproteinase 17 (ADAM17), plays a key role by cleaving membrane-bound TNFα (mTNFα) to release soluble TNFα. In normal circumstances, cellular homeostasis is maintained by ADAM17-mediated cleavage of the membrane-bound TNFR1 (mTNFR1) to release its extracellular soluble form (sTNFR1), that can bind TNFα to block inflammation [16]. However, ADAM17 has ∼100-fold higher affinity for the mTNFα compared to mTNFR1 [17]. ADAM17 is secreted by various cell types including lung cancer cells, but following RT-mediated tumor cell killing, the possibility exists for limiting levels of active ADAM17, which along with differential affinity would lead to a decrease in sTNFR1 yet will maintain higher levels of soluble TNFα to promote inflammation. Prior to the current study, the involvements of ADAM17 in TNFα-TNFR1 mediated inflammation during RP have not been studied in detail.

The involvements of multiple other cytokines, particularly IL-4, -8, -12, -13, TGF-β1, MCP-1 (CCL-2), have been described in preclinical RP models as well as patients showing pneumonitis [Ref. [18] and references therein]. In addition, studies reporting the roles of circulating miRNA involved in RP have started to emerge. In one of the earlier studies, Dinh et al. identified decreased expression of miR-29a-3p and miR-150-5p following RT in NSCLC patients [19]. In addition, the polymorphism of miR-125a, which lowers its expression, has been linked to RP [20,21].

RP is a complex disease likely to involve the interaction of various lung cell types, resident and infiltrating immune cells and tumor cells. All these cells release multiple cytokines, chemokines, miRNAs and other factors, whose expressions are further altered following RT. However numerous investigations of biological markers, predominantly using either cytokines or miRNAs, have not produced a consensus on the inflammatory pathway(s) involved. Data-driven Bayesian network (DD-BN) computer algorithms can handle multiple parameters at once, and we have previously utilized such an approach to successfully predict local tumor control in NSCLC patients for response-adopted radiotherapy [22]. In the current study have combined several approaches to analyze biological variables (circulatory cytokines and miRNAs collected prior to and multiple times during RT) and examine their relationships to mean lung dose (MLD) and clinical parameters to identify a biomarker signature for NSCLC patients who may be at higher risk of developing RP2.

Materials and methods

Clinical trial design

Samples were collected from a pilot prospective study of response-driven adaptive radiation therapy for patients with locally advanced NSCLC (NCT02492867). This study was designed to apply functional imaging, Fluorodeoxyglucose-Positron Emission Tomography (FDG-PET) ("a PET scan") and Ventilation/Perfusion Single Photon Emission Computerized Tomography (V/Q SPECT) ("a perfusion scan") before treatment and then again during treatment to study therapeutic efficacy and lung function. Blood samples were collected for biomarker studies pre-treatment, during-treatment (at day 1, 2, 5, week 2, 3, 4, 5 and 6 during chemoradiation), and in follow-up after chemoradiation (1-month post-RT, and then every 3-months post-RT for the first year). Blood samples were immediately processed in the laboratory after collection, as described below. Patient were followed clinically after concurrent chemoradiation at 1 month and 3 months post-treatment, and then every 3 months for the first year, then every 6 months until year 3, then annually through year 5. Toxicity, including pneumonitis, was recorded and graded according to Common Terminology Criteria for Adverse Events (CTCAE v4.0).

Radiation therapy

Patients received radiation treatment 5 days per week, in once daily fractions, for 30 treatments with dose per fraction individually adapted over the final 9 treatments based on individual FDG-PET response and V/Q SPECT response during treatment [23]. The intent of delivering response-based treatment was to intensify the dose to active tumor while limiting normal tissue toxicity to an estimated RP2 rate of 15%, and grade 2 or greater esophagitis rate of 30% [24]. Patients also received concurrent chemotherapy with carboplatin and paclitaxel. Following completion of concurrent chemoradiation, patients received consolidation chemotherapy (carboplatin and paclitaxel) or ICB (durvalumab) at the discretion of the medical oncologist.

Specimen collection

All specimens were prospectively collected under an institutional review board-approved protocol (HUM00098202) after receiving informed patient consent. Plasma samples were processed from blood collected in EDTA tubes (Fisher, Cat No. 02-657-32). Briefly, within one hour of blood draw, the blood was centrifuged at 3400 × g for 10 min at room temperature. The supernatant was then centrifuged at 1940 × g for 10 min at room temperature with no brake during the deceleration. The resulting supernatant was then aliquoted and frozen at −80 °C.

RNA isolation from patient plasma

Patient plasma (same time frame collected from different patients) was thawed at 37 °C for a maximum of 3 min, with a batch limit of 12 samples handled at once. RNA was extracted utilizing the miRNeasy Micro Kit (Cat#217004, Qiagen) according to manufacturer's instruction. In brief, 200 μl of plasma was mixed by vortexing with 1 mL of Qiazol, incubated at room temperature for 5 min, then mixed with 200 μl of chloroform, and centrifuged. The upper aqueous phase was mixed with 1.5 volumes of ethanol and purified using the miRNeasy mini column. Isolated miRNAs were stored at −80 °C until analysis.

Quantification miRNA using OpenArray

Isolated RNA was first reverse transcribed and pre-amplified using separate A and B primer pools, following the Life Technologies OpenArray protocol (4461306 Rev. C) with modifications for low input (Life Tech Application Note). TaqMan OpenArray Human MicroRNA, QuantStudio 12K Flex Panels (part no. 4470187) and the necessary reagents indicated in the OpenArray protocol were provided by Life Technologies. The Accufill robot was used to load panels in accord with standard protocols (part no. 4456986 Rev. A), and the sample was thermally cycled immediately in the QuantStudio 12K Flex Real-Time PCR instrument. Raw data (*.eds files) were exported in *.txt files for analysis with Microsoft Excel and other statistical R software packages.

Multiplex cytokine analysis

Levels of selected plasma cytokines from samples collected at pre-, 2 and 4 weeks (represent the most complete dataset) during-treatment time points were quantified using Procartaplex 15-plex customized ELISA kit (Cat# PPX-15, Life Technologies, Grand Island, NY) according to manufacturer's instruction. In brief, ProcartaPlex™ immunoassays use magnetic microsphere technology to detect multiple protein targets, here cytokines and chemokines. In the current study, based on our prior findings and literature summary, we selected 15 target inflammatory cytokines [sTNFR1 (for TNF-α), IFN-γ, IL-1α, IL-1β, IL-4, -6, -8, -15, -17α, -23, IP-10 (CXCL-10), MCP-1 (CCL2), PDGF-bb, CD40L, GM-CSF] to develop a customized assay panel. Prior to running the assay, plasma samples were thawed on ice and then centrifuged at 10,000 × g. The supernatants were assayed in duplicate per manufacturer protocol and analyzed on a Luminex microplate reader.

sTNFR1 ELISA

Plasma samples were assayed for sTNFR1 using the Quantikine ELISA kit (Cat#DRT100, R&D systems, Minneapolis, MN) according to the manufacturer's protocol. Briefly, each plasma sample was diluted 1:10 in the provided diluent and assayed in duplicate. The samples were incubated, on plate, with assay diluent for 2 h at room temperature and then washed 3 times with the supplied wash buffer. Human TNFR1 conjugate was added to each well followed by another 2 h of incubation. Following an additional 3 washes, the substrate and stop solutions were added, and the plate was read on a plate reader per manufacturer protocol. Cytokine concentrations were calculated against a standard curve utilizing the TNFR1 protein standard supplied with the ELISA.

ADAM17 ELISA

Total amounts of ADAM17 present in patients’ plasma were quantified using an ELISA kit (Cat # EHADAM17, Thermo Fisher) according to manufacturer's instructions. In brief, 50 μl of plasma from each patient was incubated with Anti-Human TACE Strip Plate for 2.5 h in triplicate. The samples were processed according to the protocol, maintaining the incubation time, and washing steps. Following colorimetric assay, absorbance was measured at 450 and 550 nm, and the final values was calculated by subtracting absorbance at 550 nm from 450 nm to correct for optical imperfections in the microplate.

Immunoblotting of active ADAM17 from patient plasma

Thawed plasma samples were first diluted 50-fold with ice-cold PBS, and the protein concentrations were quantified using the Bradford method according to manufacturer's protocol (Cat. B6916, Sigma). 15 ug of protein from plasma samples was loaded in a 4–12% bis-tris gel (Invitrogen, Cat# NP0321) and immunoblotted as standardized and described earlier [15]. For the detection of active ADAM17, we purchased an antibody (Cat# ABT94, Sigma) at a dilution of 1: 50,000.

Statistical analysis

miRNA analysis

To obtain the most reliable miRNA values and their associated changes, we performed three sequential steps: quality assessment, non-specific filtering and differential expression analysis as described in a recent publication [25]. In brief, a miRNA 1:4 dilution experiment was performed to determine which miRNAs had cycle threshold (Ct) value measurements consistent with the expected changes from the dilution. Of these reliable miRNAs, those that were not variant across samples were removed, as they would not be effective at predicting RP. Patients (n = 26) were divided into pneumonitis and non-pneumonitis groups, and the CT values of the remaining miRNAs were compared across groups using t-tests, to assess the individual significance of each miRNA. We used false discovery rate adjusted p-values to account for multiple hypothesis testing. We also sought to predict RP2 by considering multiple miRNAs simultaneously, which was performed through multivariable linear regression. Models were selected via a forward selection process using the Akaike Information Criterion (AIC), and variables with high variance inflation factors were removed to address multicollinearity issues. This process was applied three times to create models using pre-treatment miRNA values, during treatment miRNA values, and delta (pre-treatment minus during treatment) values as predictors for RP.

Cytokine analysis

For cytokines, we applied univariate linear, multivariate linear regression and log transformation analyzes. Based on univariate analysis, none of the cytokines showed statistical significance. Also, due to small sample size (n = 39) with a larger number of variables (15 covariates), Lasso was impossible to perform. For stepwise regression for feature selection based on AIC (Akaike Information Criterion), we first developed a model (Model 1) by fitting linear regression model with all covariates. We then conducted stepwise regression for Model 1 to get Model 2, followed by calculation of variance inflated factor (VIF) of Model 2 and deleted the variable with highest VIF. We repeated the above steps until all the VIFs were less than 10.

Data-driven Bayesian network (DD-BN) analysis

To unravel radiobiological interactions of cytokines and miRNA biomarkers along with the TNFα-NFκB signaling pathway in RP2 prediction, our previously developed graphical DD-BN approaches were employed to identify the most robust radiobiological variables related to RP2 prediction [22,26]. A BN is a probabilistic directed acyclic graphical model that uses Bayesian inference to estimate conditional dependence of involved variables. Our DD-BN approach further allows exploring these variables' potential causal-effect relationships in terms of response to radiation treatment through robust BN structure learning and parameter estimation before and during radiotherapy. As radiation induced lung toxicity is influenced by dose per fraction [27], for DD-BN modeling we decided to use gEUD [22,26] instead of mean lung dose (MLD), a more widely used parameter directly obtained from radiation therapy planning. Nested cross-validation and area under receiver operating characteristic (ROC) were used to evaluate the prediction performance of these DD-BN models. The prediction 95% confidence intervals were calculated using statistical resampling with 2000 stratified bootstrap replicates.

Results

Patient characteristics

This study included 39 locally advanced non-small cell lung cancer patients, and the characteristics (demographic clinical and treatment data) are shown in Table 1. In this cohort, we recorded 8 cases (∼20%) of grade 2 and above radiation pneumonitis cases.

Table 1

Demographic, clinical, and treatment data.

Parameters	All Patients (n = 39)
Age (years)
Median (range)	66 (54, 81)
Sex, n (%)
Male	26 (66.7%)
Female	13 (33.3%)
Zubrod Performance Status, n (%)
0	8 (20.5%)
1	29 (74.4%)
2	2 (5.1%)
Smoking Status, n (%)
Never	2 (5.1%)
Former	24 (61.5%)
Current	13 (33.4%)
Stage Group, n (%)
IIA	1 (2.6%)
IIIA	18 (46.1%)
IIIB	18 (46.1%)
IIIB-C	1 (2.6%)
IIIIA	1 (2.6%)
Concurrent Chemotherapy, n (%)
Yes	39 (100%)
Radiation Dose, Gy
Physical Dose
Median (range)	67.8 (59.97, 80.40)
EQD2
Median (range)	69.29 (60.20, 86.30)
BED (Gy)
Median (range)	83.15 (72.24, 103.56)
Mean Lung Dose (Gy)
Median (range)	15.131 (5.745, 21.1)
Highest Pneumonitis Grade n (%)
0	24 (61.6)
1	7 (17.9)
2	7 (17.9)
3	1 (2.6)

Alteration of circulating cytokines that predict RP

Mean lung dose (MLD) is a known confounder to consider while analyzing the role of cytokines and miRNA on RP [28]. To understand the proportion of variation explained by MLD, a linear regression model showed that in our dataset MLD alone was not sufficient to predict toxicity (p-value = 0.125 > 0.05). Therefore, we decided to include MLD when fitting the subsequent multivariate regression model using cytokines as a biomarker. We performed a customized multiplex ELISA screening for 15 selected inflammatory cytokines (as described in Materials and Methods) in 39 patients' plasma samples (with 8 reported toxicities). Table 2 shows the association between the lung toxicity and the cytokines via multivariate linear regression model. The cytokine levels used for analysis were either prior to radiation therapy (base line), or after 2- or 4-weeks during treatment. The model accuracy is shown by the coefficient of determination as adjusted R2. For the base line cytokine values, the adjusted R2=0.222 after leave-one-out cross-validation (Table 2A), which included 7 cytokines (IL-4, IP-10, IL-6, IFN-γ, MCP-1, IL-23, and PDGF-bb). Cytokine expression differences noted at week 2 (IL-4, -6, -1α with R2 value of 0.183) and week 4 (IP-10, IL-8 and -23 with R2 value of 0.219) (Table 2B,C). The cytokine signature obtained based on the difference between week 2 and week 4 from the baseline and between week 4 minus week 2 are included in Table 2D–F, respectively.

Table 2

Multivariate linear regression model showing identified cytokine signature those correlate with highest pneumonitis (RP) grade at the indicated time points.

A. Fitting model between pretreatment cytokine values and highest RP grade (adjusted R2 = 0.222)
Cytokines	Estimate	Std. Error	Pr (>ItI)
IL-4	0.645	0.237	0.011
IP-10	0.333	0.203	0.110
IL-6	0.422	0.187	0.031
IFNγ	−1.354	0.247	0.00001
MCP1	−0.399	0.120	0.002
IL-23	−0.483	0.167	0.007
PDGF-bb	0.685	0.252	0.011
MLD	0.066	0.026	0.019
B. Fitting model between Wk2 cytokine values and highest RP grade (adjusted R2 = 0.183)
Cytokines	Estimate	Std. Error	Pr (>ItI)
IL-4	0.430	0.275	0.130
IL-6	−0.885	0.285	0.005
IL-1α	0.493	0.216	0.031
MLD	0.094	0.037	0.018
C. Fitting model between Wk4 cytokine values and highest RP grade (adjusted R2 = 0.219)
Cytokines	Estimate	Std. Error	Pr (>ItI)
IP-10	0.435	0.123	0.001
IL-8	−0.968	0.222	0.0001
IL-23	0.694	0.209	0.002
MLD	0.086	0.030	0.008
D. Fitting model between Wk2 minus baseline cytokine values and highest RP grade (adjusted R2=0.183)
Cytokines	Estimate	Std. Error	Pr (>ItI)
IL-8	0.698	0.312	0.035
IL-23	−0.440	0.263	0.108
IL-15	0.675	0.406	0.110
CD40L	−0.707	0.266	0.014
PDGF-bb	−0.573	0.307	0.075
MLD	0.088	0.038	0.032
E. Fitting model between Wk4 minus baseline cytokine values and highest RP grade (adjusted R2=0.199)
Cytokines	Estimate	Std. Error	Pr (>ItI)
IP-10	−0.444	0.124	0.001
MCP-1	0.262	0.121	0.038
IL-1α	−0.328	0.118	0.009
MLD	0.086	0.030	0.008
F. Fitting model between Wk4 minus Wk2 cytokine values and highest RP grade (adjusted R2 = 0.3)
Cytokines	Estimate	Std. Error	Pr (>ItI)
IL-1b	0.919	0.214	0.0005
IP-10	−0.476	0.255	0.079
sTNFR1	−0.581	0.197	0.009
IFNγ	0.609	0.309	0.065
GM-CSF	−0.561	0.246	0.036
IL-1α	−0.195	0.153	0.220
IL-23	−0.575	0.208	0.013
MLD	0.118	0.030	0.001

High basal level of sTNFR1 and its transient drop correlate with RP

After identifying inflammatory cytokine signature as potentially predictive biomarkers of radiation-induced lung toxicity, we wished to assess the involvement of TNFα. We previously identified the importance of TNFα signaling in RP in a mouse model [15], and here our cytokine signature from patient samples identified a drop in sTNFR1 as one of the predictors of RP. Recognizing, as described in the introduction, that sTNFα, which is very labile, and sTNFR1 are in a dynamic balance it was of interest to quantify sTNFR1 [29] both at the base line and at multiple time points during RT. As shown in Fig. 1A, prior to radiation therapy, the average sTNFR1 level was 1334.8 ± 585.6 pg/mL, which is significantly higher than previously reported values found among healthy individuals (∼200 pg/mL) [30], suggesting the existence of basal inflammation in the majority of our lung cancer patients. While analyzing plasma samples collected from matched patients at multiple time points (day 1, 2, 5, week 2, 3, 4, 5 and 6) during RT, we noted RP symptoms in patients who showed lowering of sTNFR1 within 2–4 weeks of RT (Fig. 1B).

Fig. 1

High basal level of TNFα and a transient decrease in sTNFR-1 levels correlate with RP. (A) Soluble TNFR1 (sTNFR1) levels present in the patients’ plasma (n = 38) prior to RT was calculated using sandwich ELISA. The red (1334.8 ± 585.6 pg/ml) dotted line indicates average levels of sTNFR1 in lung cancer patients. The black dotted line indicates sTNFR1 levels among healthy individuals as reported in a prior study [30]. (B) sTNFR1 levels were calculated in patients’ plasma collected pre-treatment (baseline), and after different time points following initiation of RT (day 2, 5, week 2, 4, and 6), and after 1, 3 and 6 months of 5-FU treatment. Log ratio relative to the baseline sTNFR1 [logΔ(sTNFR1)] was plotted against day post first fraction of RT received. A logistic model of predicting the occurrence of RP2 toxicity was fit based upon logΔTNFR1(14), i.e. the estimated log-ratio of sTNFR1 comparing 14 days after the first fraction to pre-treatment. The assumption of these models is that all toxicities occurred after 2 weeks’ time. The blue and red thick lines give the group-wise (no RP2 vs RP2+, respectively) smoothed average with corresponding confidence band. The vertical line marks day 14, against which the logistic models were fit.

Change in ADAM17 levels in RP and non-RP patients

Active ADAM17 is involved in shedding of both the membrane bound TNFα and TNFR1, but with differential affinity; ADAM17 has 100-fold higher affinity for membrane-bound TNFα as compared to TNFR1 [17]. Hence, we compared the levels of ADAM17 in non-RP and RP patients with the hypothesis that in an ADAM17 low environment, sTNFR1 (anti-inflammatory) shedding may be compromised, while circulating TNFα (inflammatory) may remain high to promote inflammation noted among RP patients (Supplementary Fig. S5). Ionizing radiation is reported to activate ADAM17 [31], however, we identified two patterns for ELISA measuring ADAM17 protein levels in the plasma of NSCLC patients undergoing chemoradiation as part of this observational clinical trial. The first showed a subset of patients with initially very low levels of ADAM17, where <10 pg/mL average ELISA units were chosen as the cleanest visual threshold. These samples represented 20.5% of patients (8/39) with similar proportions in RP (3/8, 37.5%) and no-RP (5/30,16.7%) patient groups (Fisher Exact P = 0.15; Fig. 2A). In these samples ADAM17 levels tended to remain low following irradiation, with 4 samples in particular (50%) remaining below the 10 pg/mL threshold for both WK2 and WK4 measures, as opposed to 1/5 (20%) and 3/25 (12%) in RP and non-RP higher-expressing groups, respectively (Fig. 2B,C). Since the ELISA we performed cannot distinguish active and total protein forms, the ultra-low levels can be taken as low ADAM17 activity, while higher levels may not always correlate with higher activity as multiple factors are known to influence ADAM17 catalytic activity [32]. Kruskal-Wallis non-parametric comparisons (Dunn corrected for three sample groups) for WK2 and WK4 time points showed no significant differences in total ADAM17 levels between RP and non-RP groups (Fig. 2A–C).

Fig. 2

RP patients showed lowered active ADAM17 levels following RT. Total ADAM17 concentration (pg/mL) measures for 38 NSCLC patients, represented as a log2 scale for (A) before treatment (pre-treat), (B) at week2 (WK2) and (C) at week 4 (WK4) of radiation treatment, as described in Methods. Each datapoint represents the mean of 3 independent ELISA reactions on the same patient plasma sample. Samples were divided into 3 groups based on a) low initial levels (<10 pg/ml) indicated by triangular symbols vs those with higher levels (>10UI/ml) circular symbols, as well as to distinguish RP2 (red symbols) from those of patients that did not develop RP (non-RP; blue symbols). The only significant 2 group post-hoc comparisons are shown, based on an Kruskal-Wallis non-parametric tests at WK2 and WK4 timepoints, with Dunn adjustments for multiple testing. In WK4, the comparison between samples with low initial concentration, and the non-RP sample group yielded P = 0.0097. (D) Representative immunoblots of patients’ plasma collected pre-treatment, 2 and 4 weeks during treatment showing active ADAM17 levels among non-pneumonitic and pneumonitic patients. For normalization, equal amounts of protein (as quantified using Bradford method) were loaded. (E) Band intensities were quantified using ImageJ considering pre-treatment active ADAM17 level as ‘1’. Representative immunoblotting and quantification were performed for the RP (n = 6) and non-RP (n = 18) samples showing changes in active ADAM17 levels at week 2 and week 4 during RT. Welch's t-test was performed to calculate statistical significance. Analyzes of (F) WK2 and (G) WK4 fold changes (FC) relative to baseline (pre-treat) band intensities, considering 3 groups; very low baseline pre-pro-ADAM17 levels vs those with higher baselines that either did or did not develop RP. Data indicate a difference whereby RP samples with moderate to high baseline values tend to show reduced levels at WK2 (P = 0.048 vs low group and P = 0.025 vs non-RP). WK4 shows a similar, not nonsignificant trend. Meanwhile both non-RP, and samples with initially low pre-pro-ADAM17 levels show a similar propensity for upward trending values in WK2 and WK4. P values provided by Kruskal-Wallis tests with Dunn adjustment for multiple 2-group comparisons.

To focus specifically on functional ADAM17 levels, we assessed plasma levels for active ADAM17 using immunoblotting for 24 of the 39 total samples and using fold change (FC) to compare band intensities for WK2 and WK4 to those of baseline (pre-treatment) for each available sample (Fig. 2D,E). When FC was considered for the three sample groups (low ADAM17 ELISA; n = 6; as well as RP and non-RP with higher ELISA levels; n = 7 and 17, respectively), there was a trend towards lower FC levels for RP vs non-RP samples across both weeks, with a Kruskal-Wallis P = 0.025 at WK2 (Fig. 2F) and P = 0.23 at WK4 (Fig. 2G). Interestingly FC levels for both low ELISA and non-RP sample groups showed similarly high FC values, averaging over 1 (generally increased) for both WK2 and 4 (Fig. 2F,G), recapitulating prior data [31]. Collectively, the above data suggest that a subset of samples have consistently low and less amount of active ADAM17, while a subset with higher total ADAM17 levels with limited catalytic activity in the early stages of radiation exposure.

Alteration of circulating miRNA that predict RP

To strengthen biomarker signature, besides cytokines we also quantified plasma circulating miRNA. To determine the assay sensitivity and to identify most reliable miRNAs, we initiated our studies by quantifying equimolar and ratiometric pools of chemically synthesized RNA oligonucleotides as described previously [25]. The raw data quantifying miRNA consisted of CT values for 818 miRNAs from 26 out of total 39 patients (2 RNA samples failed quality control and 11 others were collected post analysis) taken prior to RT (baseline) and at WK4 during RT. The list of 70 reliable miRNA (Supplementary Table 1) and a boxplot of these un-normalized CT values is shown in Supplementary Fig. S1. In the dilution experiment, we analyzed RNA isolated from 3 patients and used two replicates of 1:4 dilutions. For a reliable miRNA, we expected the difference in CT value from the quarter dilution sample to the undiluted sample to be 2. We accepted a miRNA as reliable if, for each run of the experiments (3 patients, 2 replicates each), the difference in CT value for the that miRNA was between 1 and 3. This narrowed down the initial pool of 818 miRNAs to 70. A boxplot consisting of the 70 reliable miRNAs is shown in the Supplementary Fig. S2. For each miRNA, we calculated the variance of the delta (pre-treatment minus during treatment) values across all the patients and removed the miRNAs in the bottom 10% of delta value variances. This filtering reduced the number of variables to 63 miRNAs, of which 9 were U6 miRNAs; thus, we had 54 miRNAs to analyze.

For each of the 26 patients, we had an outcome of highest pneumonitis grade, denoted as 0, 1, 2, and 3, with 3 being the most severe toxicity. To create pneumonitis and non-pneumonitis groups, patients with the highest grade of 0 and 1 formed the non-pneumonitis group, while 2 and 3 formed the pneumonitis group. The histogram of the p-values of the delta value t-test showed potential signs of significance (Supplementary Fig. S3A). However, after performing false discovery rate calculations (q-value), all the miRNAs failed to be significant (Supplementary Fig. S3B). This suggested that the differences in individual miRNA CT values between the pneumonitis and the non-pneumonitis groups were not significant.

We then created models involving linear regression and stepwise feature selection with highest pneumonitis grade as the outcome, and the 54 miRNAs as predictors. For the forward selection approach, we used the Akaike Information Criterion (AIC) to select miRNAs to predict pneumonitis. For the pre-treatment miRNA value model, we obtained 5 significant miRNAs (miR-532, -451, -99b, -139-3p, and Let-7c) with an adjusted R2 of 0.6385 (Table 3A). For the delta value model (WK4 minus pre-treatment), we obtained one significant miRNA (miR-495) with an adjusted R2 of 0.3386 (Table 3B). Since the predictors are the CT values, and CT values are inversely proportional to the amount of miRNA, a positive coefficient represented a negative association between the miRNA and toxicity, and vice versa. From the above analysis we identified three miRNAs (miR-532, -495, -99b) positively correlated with RP2, whereas the remaining three (miR-139-3p, -451, and let-7c) were negatively correlated. Interestingly, reported studies on these aforementioned miRNAs identified them as regulators of NFκB [33], [34], [35], [36], [37], [38], a key mediator of multiple inflammatory signaling processes, including radiation-induced normal tissue toxicity [39].

Table 3

Linear regression and stepwise selection with highest Pneumonitis grade identified selected miRNAs pre-treatment (in A) and Wk4 during treatment (in B) signature using Akaike Information Criterion (AIC).

A. Fitting model between pre-treatment miRNA values and highest RP grade (R2 = 0.6385)
miRNA	Estimate	Std. Error	Pr(>ItI)
miR-532	−0.133	0.040	0.003
miR-451	0.316	0.075	0.0004
miR-99b	−0.405	0.083	0.00009
Let-7c	0.222	0.059	0.001
miR-139-3p	0.099	0.055	0,08
B. Fitting model between pre-treatment miRNA values and highest RP grade (R2 = 0.3386)
miRNA	Estimate	Std. Error	Pr(>ItI)
miR-495	−0.377	0.086	0.0002

A data-driven Bayesian network (DD-BN) approach to predict RP

We have previously developed a statistical method based on resampling and information theory to guide robust development of BN structure learning, where we integrated radiobiological variables in a statistical resampling framework using bootstrap generated samples of the data to represent their interaction graphically and allow further inferences [22,26]. Utilizing the candidate variables from the statistical analyzes, we generated stable BN structures using either miRNA, cytokines or both and combined them with radiobiological and clinical parameters to best represent RP2 prediction, as shown in Fig. 3A–F and Supplementary Fig. S4. In these BNs, arcs with green or red color represent positive or negative association between any two variables in the graph, while arcs with gray color indicate a mixture of positive and negative association depending on the expression state of the involved variables. The area under the ROC curves (AUCs) show that the performance of these BNs for RP2 prediction improved significantly when all parameters were combined achieving an AUC of 0.87 (95% 0.75–0.97) on cross-validation testing.

Fig. 3

Pre- and during-treatment DD-BN for improved prediction of RP2. (A–C) Left panels represent pre-treatment BN models using cytokines (in A), miRNAs (in B) and both (in C) along with dosimetry parameter (lung_gEUD). Right panels show corresponding ROC curves of pre-treatment DD-BN based on internal cross-validation. (D-F) Left panels represent 4 weeks during-treatment BN models as above using cytokines, miRNAs and both respectively and right panels show respective ROC curves similarly analyzed and validated as above.

Discussion

In this study, we have adopted an approach using data-driven Bayesian Networks to understand the multiple interactions of miRNA array and multiplex cytokines analysis, along with radiobiological (such as MLD and gEUD) and clinical data to predict RP2. We identified specific miRNAs (miR-532, -451, -99b, -139-3p, -495 and let-7c) and found a signature using the WK4 minus WK2 difference which included 7 cytokines (IL-1β, IP-10, IL-23, sTNFR1, IFN-γ, GM-CSF, and IL-1α) (Fig. 2B). This signature included either regulators (TNFα and IL-1α) or targets of NFκB, a key downstream mediator of TNFα-driven inflammation [39]. We further identified higher basal levels of sTNFR1 in majority of our lung cancer patients, however, surprisingly we observed a temporal (2-4 weeks during RT) reduction of anti-inflammatory sTNFR1 in RP2 patients, suggestive of a persistent inflammatory state. Using DD-BN analysis, we combined biological, radiobiological, and clinical data, which highlight the interactions between these variables as discussed below.

A key finding of our work is that it required a combined assessment of miRNAs, cytokines along with MLD to best predict radiation induced lung inflammation indicating crosstalk between multiple biomarkers. We previously reported that lung irradiation induces early release of TNFα, and inhibition of TNFR1-mediated inflammatory signaling can radioprotect mouse lung [14]. Conversely, Tristetraprolin (TTP, Zfp36 in mice) knockout mice with high basal levels of TNFα were prone to RP [14,15,40]. Such preclinical data allowed us to predict conservation of the TNFα-TNFR1 pathway in clinical settings. Analysis of patients' plasma using a multiplex cytokine array (Table 2) and conventional sandwich ELISA (Fig. 1), showed that transient loss of sTNFR1 is one of the contributing factors promoting RP. Besides activation of TNFα signaling, involvements of IL-1α and IL-6 have been reported in independent clinical studies [41,42], which is consistent with our findings. In another study with NSCLC patients during and after radiotherapy, changes in the levels of IP-10, MCP-1 and IL-6 were found to be correlated with MLD [43], similar trend was noted in our study. However, to better predict patients prone to develop RP, multiple biomarker analyzes may be required where our approach of using DD-BN analysis will be critical.

A possible involvement of TNFα-NFκB signaling in RP became evident while analyzing miRNA data. From unbiased screening of the total 818 miRNAs, we identified six [6] (miR-532, -99b and -495, let-7c, miR139-3p and miR-451), which were correlated with RP incidence. Interestingly, these miRNAs were previously reported to regulate NFκB signaling, either as activators or inhibitors. Among activators, miR-99b has been reported to promote NFκB-mediated inflammation via activating expression of IL-1β, -6 and -12 [35]. Similarly, miR-495 is known to be involved in inflammatory diseases, particularly Crohn's disease, where it regulates the expression of the nucleotide binding oligomerization domain containing protein 2 (NOD2), a key inducer of multiple cytokines including IL-6, -8 and TNFα (34). For expression of these inflammatory cytokines NOD2 is an important activator of NFκB signaling [44]. Although the direct role of miR-532 is not known, its connection with NFκB signaling results via regulation of KIFC1 expression to impact gankyrin/AKT signaling [45] and thus can indirectly impact NFκB; at the same time, the promoter of miR-532 has an NFκB1 binding site for direct regulatory feedback [33]. This miRNA is also aberrantly expressed in patients with inflammatory bowel disease [46]. In contrast, miR-451, -139-3p and let-7c are involved in NFκB regulation by negatively impacting transcription activation. For example, miR-451 directly suppresses IKKβ expression, the critical kinase that inactivates the inhibitor of NFκB signaling, IκB, via phosphorylation to promote proteasomal degradation [47]. This miRNA also suppresses TLR-4-mediated expression of inflammatory cytokines including IL-6, -1β and TNFα. Downregulation of this miRNA is also reported in nonalcoholic steatohepatitis [48]. Similarly, miR-139 is known to inhibit NFκB-mediated expression of IL-6, 1β and TNFα; thus, its targeted deletion promoted intestinal inflammation and colorectal cancer in a mouse model [36]. Furthermore, let-7c expression is reported to be lower in chronic obstructive pulmonary disease (COPD) patients [49], and its anti-inflammatory role in LPS-induced inflammation is mediated via suppression of STAT3 expression, a known collaborator of NFκB [50]. Taken together, the current literature supports our miRNA and cytokine data, which identified a signature in lung cancer patients undergoing thoracic radiation that resulted in the activation of TNFα-NFκB signaling to cause RP.

One of the major strengths of this study is the use of DD-BN modeling, which combined biological data (miRNAs and cytokines) with radiobiological parameters, such as lung generalized equivalent uniform dose (gEUD) [51]. As shown in Fig. 3, such an approach identified specific cytokines and miRNAs that showed either positive, negative or mixed correlation and strengthened our RP2 prediction ability both at the pre-treatment and during treatment (WK4) time points as evidenced by the improvement of AUC values. Such data further identified that besides increasing lung dose and loss of sTNFR1, increases in IL-4 and MCP-1, at WK4 of treatment, positively correlated with subsequent RP2 occurrence. A previous study using a rat model for RP showed increase in IL-4 following a single dose of 20 Gy thoracic irradiation [52]. Furthermore, a recent study identified association of IL-4 genetic variation (single nucleotide polymorphism, SNP) with prevalence of RP in lung cancer patients [53]. Similarly supportive, MCP-1 (also called CCR2) is known to be altered in multiple inflammatory diseases, particularly following thoracic irradiation [54]. Thus, the DD-BN analysis helped us better understand the cross-talks between different biological mediators, and their association with lung dose to predict RP2. Such findings may have both predictive and therapeutic value, as they suggest testing of anti-TNFα agents to prevent RP2.

ICB therapy following RT has become the standard-of-care for lung cancer patients, yet the use of high dose corticosteroids to treat RP may be counterproductive due to their leukotoxicity and known negative impacts on ICB efficacy. In a recent preclinical study, early blockade of NFκB signaling using Thalidomide was shown to mitigate chronic inflammation in radiation-induced urinary bladder dysfunction [55]. Our data suggest that an anti-TNF/NFκB agent may be an ideal replacement. Thus, our signature biomarker will not only be useful in identifying patients prone to RP, but it also suggests a new therapeutic strategy to treat RP2. Such a strategy is further strengthened by recent reports and our unpublished data that demonstrate prophylactic TNFα blockade may potentiate ICB efficacy at least in preclinical models [56]. Although a recent retrospective study of patient outcomes from a Dutch National Registry identified worse outcomes in patients whose ICB toxicity was managed using anti-TNFα agent, particularly Infliximab [57], in an associated commentary, certain shortfalls were also noted with respect to data interpretation [58]. Thus, we are cautiously optimistic and believe the prevalence of data suggest that a detailed preclinical study using an anti-TNFα agent is warranted.

Another interesting finding from our study is the indication of a transient reduction of sTNFR1 level between 2-4 weeks associated with subsequent RT. sTNFR1 is the protease-cleaved extracellular domain of TNFR1, which has the ability to bind and neutralize TNFα, thus acting as an anti-inflammatory factor [59]. TNFα is initially synthesized as a membrane-bound inactive form (mTNFα), which, following an insult such as RT, undergoes processing by a metalloproteinase, called TNFα converting enzyme (TACE, also cells ADAM17) [17]. ADAM17 is also the protease for TNFR1; however, its enzymatic affinity is 100-fold higher for TNFα as compared to TNFR1 [17]. Thus, in an ADAM-17 limited state, an inflammatory signal may persist due to increased presence of soluble TNFα compared to sTNFR1. ADAM17 is initially produced in its high molecular weight inactive pro-form, which can be converted to an active-form [60]. Using immunoblotting, we noted reduced levels of active ADAM17 in patients with RP2 between 2-4 weeks during RT, suggesting a time-sensitive limit in functional ADAM17. It is interesting to note that we also found sTNFR1 shows a similar drop in RP samples at WK2. These data provide preliminary mechanistic insight (Supplementary Fig. S5), requiring more detailed temporal sample collection and analysis.

We recognize that this study has some limitations. Most importantly, the modest sample size (n = 39) and lack of validation cohort somewhat restrict our inferences. Another potential weakness is overfitting of our data, but we used statistical rigor to avoid this. Additionally, we recognize that the direct quantification of TNFα poses a persistent challenge to the field, which we have approached indirectly by measuring sTNFR1 [29] and ADAM17.

In summary, we have utilized both miRNA and cytokines and using DD-BN modeling combined biological data with radiobiological parameters to identify signature biomarkers predictive of RP2, which emphasized the role of TNFα-NFκB signaling in causing lung toxicity. As both TNFα and NFκB are known mediators of various inflammatory diseases, there are already multiple agents and FDA approved drugs targeting these molecules, which can be repurposed for possible lung radioprotection.

Data availability

All miRNA and multiplexed cytokine data generated for this study are included under Supplementary Table and Figures. The raw data will be made available by the authors, without reservation to anyone requesting.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.